Binary-black-hole initial data with nearly-extremal spins 
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There is a significant possibility that astrophysical black holes with nearly-extremal spins exist. 
Numerical simulations of such systems require suitable initial data. In this paper, we examine three 
methods of constructing binary-black-hole initial data, focusing on their ability to generate black 
holes with nearly-extremal spins: (i) Bowen-York initial data, including standard puncture data 
(based on conformal flatness and Bowen-York extrinsic curvature), (ii) standard quasi-equilibrium 
initial data (based on the extended-conformal-thin-sandwich equations, conformal flatness, and max- 
imal slicing), and (iii) quasi-equilibrium data based on the superposition of Kerr-Schild metrics. We 
find that the two conformally-flat methods (i) and (ii) perform similarly, with spins up to about 
0.99 obtainable at the initial time. However, in an evolution, we expect the spin to quickly relax 
to a significantly smaller value around 0.93 as the initial geometry relaxes. For quasi-equilibrium 
superposed Kerr-Schild (SKS) data [method (iii)], we construct initial data with initial spins as 
large as 0.9997. We evolve SKS data sets with spins of 0.93 and 0.97 and find that the spin drops 
by only a few parts in 10 4 during the initial relaxation; therefore, we expect that SKS initial data 
will allow evolutions of binary black holes with relaxed spins above 0.99. Along the way to these 
conclusions, we also present several secondary results: the power-law coefficients with which the 
spin of puncture initial data approaches its maximal possible value; approximate analytic solutions 
for large spin puncture data; embedding diagrams for single spinning black holes in methods (i) 
and (ii); non-unique solutions for method (ii). All of the initial data sets that we construct contain 
sub-extremal black holes, and when we are able to push the spin of the excision boundary surface 
into the super-extremal regime, the excision surface is always enclosed by a second, sub-extremal 
apparent horizon. The quasilocal spin is measured by using approximate rotational Killing vectors, 
and the spin is also inferred from the extrema of the intrinsic scalar curvature of the apparent hori- 
zon. Both approaches are found to give consistent results, with the approximate-Killing-vector spin 
showing least variation during the initial relaxation. 

PACS numbers: 04.25.D-,04.25.dg,04.20.Ex,02.70.Hm 



I. INTRODUCTION 

There is a significant possibility that black holes with 
nearly-extremal spins exist; by "nearly-extremal", we 
mean that the spin 5* and mass M of the hole satisfy 
0.95 < S/M 2 < 1. Some models of black-hole accre- 
tion [1~3] predict that most black holes will have nearly- 
extremal spins, and observational evidence for black holes 
with nearly-extremal spins includes, e.g., estimates of 
black-hole spins in quasars [4] and estimates of the spin 
of a black hole in a certain binary X-ray source [5] . There 
is considerable uncertainty about whether black holes do 
in fact typically have nearly-extremal spins; e.g., some 
models [6-8] of black-hole accretion do not lead to large 
spins. This uncertainty could be reduced by measuring 
the holes' spins directly using gravitational waves. 

This prospect of detecting the gravitational waves 
emitted by colliding black holes, possibly with nearly- 
extremal spins, motivates the goal of simulating these 
spacetimes numerically. Indeed, one focus of intense re- 
search has been spinning black hole binaries, including 
the discovery of dramatic kicks when two spinning black 
holes merge [9-17] as well as some initial exploration of 
the orbital dynamics of spinning binaries [18-23]. All 
of these simulations start from puncture initial data as 
introduced by Brandt and Brugmamr [24] . 



The simplifying assumptions employed in puncture ini- 
tial data make it impossible to construct black holes with 
spins arbitrarily close to unity. The numerical value of 
the fastest obtainable spin depends on which dimension- 
less ratio is chosen to characterize "black hole spin." Of- 
ten, dimcnsionless spin is defined based on quasilocal 
properties of the black hole, 



where S is taken to be nonncgativc and is a suitable 
quasilocal spin (e.g., obtained using approximate rota- 
tional Killing vectors on the apparent horizon as de- 
scribed, for example, in Appendix A) and M is a suit- 
able quasilocal mass. The latter may be obtained from 
Christodoulou's formula relating spin, area and mass of 
a Kerr black hole, 

M2: = M i« + zSr> ( 2 ) 

irr 

where we define the irreducible mass in terms of the area 
A of the apparent horizon by M- m := A/IQtt. 

The quantity \ is not preserved during an evolution. 
Specifically, most black hole initial data are not exactly 
in equilibrium, which leads to transients and emission 
of an artificial pulse of gravitational radiation early in 
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numerical simulations. The geometry in the vicinity of 
the black holes relaxes on a time-scale i rc i ax (typically a 
few M), and during this relaxation, the spin changes by 

Ax := X (t = 0) - X (irdax) ■ (3) 

When constructing a single spinning black hole with 
standard puncture data [24], for instance, %(i = 0) < 
0.98, which seems encouragingly large. However Dain et 
al. [25, 26] evolved standard puncture data with initial 
spin close to this limit, and they find that the spin rapidly 
drops to x(*reiax) ~ 0.93, i.e. Ax ~ 0.05. 

For single-black-hole spacetimes, another widely used 
dimensionless spin-measure is the ratio of total angular 
momentum 1 Jadm and Arnowitt-Deser-Misner (ADM) 
energy £ A dm , 



Dain et al. noted that x(^reiax) is close to sj and ex- 
plained this result as follows: the spacetime is axisym- 
metric, which implies that the angular momentum Jadm 
is conserved and that the black hole's spin equals Jadm- 
Moreover, so long as a negligible fraction of the space- 
time's energy is carried off by the spurious radiation, the 
hole's quasi- local mass will relax to a value of -Eadm, giv- 
ing x (trciax) ~ £J- Thus conformally-flat Bowcn-York 
data cannot be used to simulate black holes with nearly- 
extremal equilibrium spins, even though the initial spins 
can be made fairly close to x = 1- 

This paper examines three different approaches of con- 
structing black hole initial data with nearly-extremal 
spin. First, we revisit puncture initial data and inversion- 
symmetric Bowen-York initial data. We show that for a 
single, spinning black hole at rest, both approaches are 
identical, and we determine spin-limits based purely on 
initial data more accurately than before: 

e.j < 0.928200, X (t = 0) < 0.9837. (5) 

We show that the limiting values of ej and x(i = 0) 
are approached as power-laws of the spin-parameter (cu- 
riously, with different powers). We furthermore give 
insight into the geometric structure of these high-spin 
Bowen-York initial data sets through numerical study 
and approximate analytical solutions and find that a 
cylindrical throat forms which lengthens logarithmically 
with the spin-parameter. 

Second, we investigate the high-spin limit of another 
popular approach of constructing initial data, the quasi- 
cquilibrium formalism [27-31] based on the conformal 
thin sandwich equations [32, 33]. For the standard 
choices of conformal flatness and maximal slicing, we are 



1 We define here Jadm by an ADM-like surface integral at infin- 
ity; in axisymmetry this definition coincides with the standard 
Komar integral for angular momentum (see Sec. II B for details.) 



able to construct initial data with spins somewhat larger 
than the standard Bowen-York limits given in Eq. (5): 

sj < 0.94, x (t = 0) < 0.99. (6) 

Once again ej is much lower than \{t — 0), which sug- 
gests that these data sets lead to equilibrium spins of 
approximate magnitude x ~ 0.94. Interestingly, these 
families of initial data are found to exhibit non-unique 
solutions [34-36] , and the largest spins are obtained along 
the upper branch. 

The third approach also utilizes the quasi-equilibrium 
formalism [27-31], but this time we make use of the free- 
dom to chose an arbitrary background data. Specifically, 
we choose background data as a superposition of two 
Kerr-Schild metrics. This approach is based on the orig- 
inal proposal of Matzner and collaborators [37, 38] and 
was first carried over into the conformal thin sandwich 
equations in Ref. [39]; also, background data consisting 
of a single, non-spinning Kerr-Schild black hole was used 
to construct initial data for a black-holc-neutron-star bi- 
nary in Ref. [40]. For single black holes, this data simply 
reduces to the analytical Kerr solution. For binary black 
holes, we construct initial data with spins as large as 

X (t = 0) = 0.9997. (7) 

We also present evolutions, demonstrating that our 
rapidly-spinning initial data sets remain rapidly-spinning 
after the numerical evolution relaxes. In particular, we 
evolve an orbiting binary with X (t = 0) = 0.9275 and a 
head-on merger with x(t = 0) = 0.9701. In both cases, 
\^x/x(t = 0)| is significantly smaller than 10~ 3 . We con- 
clude that the conformally-curved SKS initial data we 
present in this paper, in contrast with conformally-flat 
Bowen-York data, is suitable for simulating binary black 
holes with nearly-extremal spins. 

Throughout the paper, we use two different techniques 
to measure the dimensionless spin of black holes, which 
are described in the appendices. The first (Appendix A) 
technique uses the standard surface-integral based on 
an approximate rotational Killing vector of the appar- 
ent horizon. We compute the approximate Killing vector 
with a variation of the technique introduced by Cook 
and Whiting [41], extended with new normalization con- 
ditions of the approximate Killing vector, and we denote 
the resulting spin "AKV spin", xakv- The second ap- 
proach (Appendix B) is based on the shape of the horizon 
in the form of its scalar curvature; specifically, the spin 
magnitudes are inferred from the minimum and maxi- 
mum of the intrinsic Ricci scalar curvature of the hori- 
zon. We call the spin inferred in this way the "scalar 
curvature spin," and we label the spin magnitudes in- 
ferred from the scalar curvature minimum and maximum 
as x'sc" and XscfS respectively. Typically, binary-black- 
hole initial data produces holes that are initially not in 
equilibrium. Therefore, we use only the AKV spin to 
measure the initial black hole spin (Sees. III-IV.) We 
use both the AKV and the scalar-curvature spin when 
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we measure the spin after the holes have relaxed to equi- 
librium (Sec. V). 

We also monitor whether any of the constructed initial 
data sets have super-extremal spins, as this may shed 
light, for example, on the cosmic censorship conjecture. 
When using the Christodoulou formula [Eq. (2)] to define 
M, the quasilocal dimcnsionless spin x is by definition 
bounded [42], x < 1- This can be seen most easily by 
introducing the parameter defined as 



equations 



c 

and then rewriting x as 



X = 1 



2Md 



(i-C) 2 
l + C 2 



(8) 



(9) 



The ratio x is therefore not useful to diagnose super- 
extremal black holes. A more suitable diagnostic is found 
in the parameter (. For Kerr black holes, the first term on 
the right-hand-side of Eq. (2) is always smaller or equal 
to the second, with equality only for extremal spin; i.e., 
£ < 1, with equality for extremal spin. This motivates 
an alternative definition of extremality [42] : a black hole 
is said to be superextremal if the second term in Eq. (2) 
is larger than the first one, i.e. if £ > 1. In this paper, we 
monitor £, which we call the spin-extrcmality parameter, 
along with the dimcnsionless spin x- We find instances 
where £ exceeds unity. Before this happens, however, a 
larger, subextremal (£ < 1) apparent horizon appears, 
enclosing the smaller, superextremal horizon (Sec. IV B, 
Fig. 12). 

This paper is organized as follows. Section II summa- 
rizes the various formalisms that we use to construct ini- 
tial data. Section III investigates single black hole initial 
data, followed by the construction of binary-black-holc 
initial data in Sec. IV. Section V presents binary black 
hole evolutions that show the good properties of super- 
posed Kerr-Schild data, and the various spin-diagnostics. 
We summarize and discuss our results in Sec. VI. Finally, 
Appendix A and Appendix B present our techniques to 
define black hole spin. 



II. INITIAL DATA FORMALISM 

Before constructing initial data for rapidly-spinning 
single (Sec. Ill) and binary (Sec. IV) black holes, we first 
summarize the initial data formalisms we will use. After 
laying some general groundwork in Sec. II A, we describe 
Bowen-York initial data (including puncture initial data) 
in Sec. II B and quasi-equilibrium cxtendcd-conformal- 
thin-sandwich data in Sec. II C. 



R 



K 2 - K i:i KV = 0, 



0. 



V, (K ij - <j"K) 



(10) 
(11) 



Here, gij is the induced metric of the slice E, with covari- 
ant derivative V,;, R := g^Rij denotes the trace of the 
Ricci-tensor Rij , and denotes the extrinsic curvature 
of the slice S as embedded into the space-time manifold 
M. 

The constraint equations (10) and (11) can be trans- 
formed into elliptic partial differential equations using 
a conformal transformation, e.g. [33]. One introduces a 
conformal metric, gij via 



9ij = '''.'A; • 



(12) 



with the strictly positive conformal factor ip > 0. Substi- 
tuting Eq. (12) into Eq. (10) yields an elliptic equation 
for ip. One furthermore decomposes the extrinsic curva- 
ture into trace and tracefree part, 



K ij =A ij + -g ij K, 
3 y 



(13) 



and splits off a longitudinal part from the tracefree ex- 
trinsic curvature, 



a 13 = —{hvy 

a 



(14) 



In Eq. (14), a is a strictly positive weight-function, the 
longitudinal operator is defined as (LV) y = 2V^V J ' — 
§<? y V kV k > and M y is symmetric and trace-free 2 . Fi- 
nally, one introduces the conformally scaled quantities 
(7 = il. ]& a 1 M % i = ifj~ 10 M l: > , which allows the momentum 
constraint [Eq. (11)] to be rewritten completely in terms 
of conformal quantities: 



A 11 = </>~ iU A y , 

A^ = -(Lvf +M ij . 
a 



(15) 
(16) 



The Hamiltonian and momentum constraints then be- 
come 



8 E ~12 J 
l„~.~ij\ 2 



~^R- ^K 2 i/j 5 + ?-A l3 A l ^~ 7 = 0, (17) 



a 



Vj -(LV) - -ip b VK + VjAf* J = 0. (18) 



Given choices for M lJ , K, g^ and <r, and also boundary 
conditions, one can solve Eqs. (17) and (18) for ijj and 
V 1 , and then assemble the (constraint-satisfying) initial 
data g^ and K l K 



A. Extrinsic curvature decomposition 



Initial data sets for Einstein's equations are given on 
a spatial hypersurface S and must satisfy the constraint 



2 It is also possible, but not necessary, to require that M IJ is di- 
vergence free. 
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Many important approaches to construct binary black 
hole initial data can be cast in this form. The various 
approaches differ in the choices for the freely specifiable 
parts and the boundary conditions. Some choices of free 
data aim for simplicity, such as Bowen-York initial data. 
Other approaches aim to preserve freedom, resulting in 
more complicated sets of equations but also more flexibil- 
ity to control properties of the resulting initial data. The 
quasi-cquilibrium cxtendcd-conformal-thin-sandwich ap- 
proach falls into this second category, and we will exploit 
precisely its inherent freedom in choosing the free data 
to construct black holes with nearly-extremal spins. 



B. Bowen-York initial data 

In this section, we describe two approaches of con- 
structing initial data based on the well-known Bowen- 
York extrinsic curvature. These two approaches, punc- 
ture data and inversion-symmetric data, differ in how 
they treat the coordinate singularity at r = 0; both 
can be obtained from the general procedure outlined in 
Sec. II A by setting a = 1, K = 0, = and by using 
a conformally flat metric 



Sij — fi 



(19) 



The momentum constraint [Eq. (18)] then reduces to 
Vj(LV) u ' = 0, which is solved by choosing the analyt- 
ical Bowen-York solutions [43, 44] . 

The Bowen-York solutions can be written down most 
conveniently in Cartesian coordinates, = 5,-,-: 



V£ = ~ [7P l + n l P k n k ] , 



Ar 
1 

~2 



—e i m b n 



(20) 
(21) 



where r = (x^x* S^) 1 ' 2 is the coordinate distance to the 
origin and n l = x l jr is the coordinate unit vector point- 
ing from the origin to the point under consideration. The 
spatially-constant vectors P l and S l parametrize the so- 
lutions 3 



Ap ~ 2r 2 



2P<V> - (S ij - n l n j ) P k n h 



(22) 
(23) 



The conformal factor ip is then determined by the 
Hamiltonian constraint [Eq. (17)], which simplifies to 



0. 



(24) 



3 In Cartesian coordinates, upper and lower indices are equivalent, 
so index positioning in Eqs. (20)-(23) is unimportant. To find 
A„ , c in another coordinate system, first compute the Cartesian 

f I 5 

components Eqs. (20)-(23), and then apply the desired coordi- 
nate transformation. 



We would like to recover an asymptotically flat space; 
this implies the boundary condition tp — > 1 as r — * oo. 

This boundary condition makes it possible to evalu- 
ate the linear ADM-momentum and ADM-like angular 
momentum of Bowen-York initial data without solving 
Eq. (24). These quantities are defined by surface inte- 
grals at infinity, 



Ji 



1 

8tt 



(K, //,, /\ ) C •-"*' dA, 



(25) 



where s l is the outward-pointing unit-normal to the in- 
tegration sphere 4 . By letting ip — > 1 in Eq. (15), one can 
replace Kij by Ay and then evaluate the resulting inte- 
grals. The choice of vector £ l determines which quantity 
is computed: For instance, £ = e x corresponds to the 
x-component of the linear ADM-momentum, £ = 3$ = 
—xe y + ye x yields the z-component of the ADM-like an- 
gular momentum 5 . For Eqs. (22) and (23), the results 



are P, 



ADM 



P l and J 



ADM 



S l , respectively. 



The ADM energy is given by the expression 

Eaum = ^ j v, (gs - 6i j g) s l dA, 



(26) 



where C/y := .</,, ./',,. Q 
Eq. (26) reduces to 



Oijg 1 -'- For conformal flatness, 



E A 



DM 



-— f d r ip dA. 
2n 



(27) 



The derivative of the conformal factor is known only after 
Eq. (24) is solved; therefore, in contrast with the linear 
and angular momenta, £adm can be computed only after 
solving the Hamiltonian constraint. 

We now turn our attention to inner boundary condi- 
tions. A l p and A l g are singular at r = 0. This singularity 
is interpreted as a second asymptotically fiat universe; 
when solving Eq. (24), this can be incorporated in two 
ways: 

• Inversion Symmetry: The demand that the so- 
lution be symmetric under inversion at a sphere 
with radius i?i nv centered on the origin [44] results 



4 At infinity, the normal to the sphere s 1 is identical to the coor- 
dinate radial unit vector n'. 

5 As is common in the numerical relativity community, we in- 
troduce the phrase "ADM angular momentum" to refer to an 
angular momentum defined at spatial infinity in the manner of 
the other conserved ADM quantities of asymptotically flat space- 
times [45], despite the fact that (at least to our knowledge), no 
such quantity is widely agreed to rigorously exist in general, due 
to the supertranslation ambiguity that exists in four spacetime 
dimensions. For recent research on this issue see [46] and refer- 
ences therein. In the present paper, this subtlety can be ignored, 
because we only compute this quantity in truly axisymmctric 
spacetimes, with £ the global axisymmetry generator, so that 
Jadm coincides with the standard Komar integral for angular 
momentum. 
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in a boundary condition for ip at r — R ln v, namely 
dip/dr = — ip/(2R lm ). The Hamiltonian constraint 
Eq. (24) is solved only in the exterior of the sphere, 
f > Rinv, and the solution in the interior can be 
recovered from inversion symmetry [44], e.g. 



Rh 



' R 2 

l inv 

-if) I — —x 



(28) 



• Puncture data: One demands [24] the appropri- 
ate singular behavior of tp for r — > to ensure that 
the second asymptotically flat end is indeed flat. 
That is, ip must behave as 



2r 



+ 1 + u(x l ) 



(29) 



for some positive parameter m p (the "puncture 
mass") and function u(x l ) that is finite and con- 
tinuous in K 3 and approaches as r — ► oo. Equa- 
tion (24) then implies an equation for u that is finite 
everywhere and can be solved without any inner 
boundaries: 



Yu 



1 



+ -f - +ur) 



(30) 



The majority of binary black hole simulations use 
puncture data, see, e.g., Refs. [9-23]. 

Both approaches allow specification of multiple black- 
holes at different locations, each with different spin and 
momentum parameters S l and P l . For puncture data 
this is almost trivial; this accounts for the popularity of 
puncture data as initial data for black hole simulations. 
In contrast, for inversion-symmetric data, one needs to 
employ a rather cumbersome imaging procedure 6 (see 
e.g. [47] for details). 

For a single spinning black hole at the origin, the ex- 
trinsic curvature A l g given by Eq. (23) is identical for 
inversion-symmetric and puncture data. For inversion- 
symmetric data, the conformal factor has the usual falloff 
at large radii, 

ij(x l ) = l + ^^- + 0(r- 2 ), asr^oo. (31) 
2r 

Using Eq. (28) we find the behavior of ip as r - * 0: 



R\, 



E A 



DM 



2i?i, 



O(r) 



as r — > 0. 



(32) 



Comparison with Eq. (29) shows that this is precisely 
the desired behavior for puncture data, if one identi- 
fies i?; nv = TOp/2 and E / {2R- lm ) = 1 + u(0). Because 



Even for a sing le black hole with P k ^ 0, Eq. (22) has to be 
augmented by additional terms of 0(r~ 4 ) to preserve inversion 
symmetry [44]. 



puncture data has a unique solution, it follows that for 
single spinning black holes, puncture data and inversion- 
symmetric data are identical, provided m p = 2Ri nv . 

For inversion-symmetric initial data for a single, spin- 
ning black hole, it is well-known [48] that the apparent 
horizon coincides with the inversion sphere, tah = -Rinv 
Therefore, we conclude that for puncture data for a sin- 
gle, spinning black hole, the apparent horizon is an exact 
coordinate sphere with radius tah = Tn p /2, despite Ag 
and u{x l ) not being spherically symmetric. 



C. Quasi-equilibrium 
extended-conformal-thin-sandwich initial data 

Another popular approach of constructing binary- 
black-hole initial data is the quasi-equilibrium extended- 
conformal-thin-sandwich (QE-XCTS) formalism [27-31]. 
Instead of emphasizing the extrinsic curvature, the con- 
formal thin sandwich formalism [32] emphasizes the spa- 
tial metric gij and its time-derivative. Nevertheless, it is 
equivalent [33] to the extrinsic curvature decomposition 
outlined in Sec. II A. The vector V 1 is identified with the 
shift /3\ 



V 



(33) 



and the weight-functions a and a are identified (up to a 
factor 2) with the lapse and the conformal lapse, respec- 
tively, 



a = 2a, 



a = 2a. 



(34) 



The tensor Ma is related to the time-derivative of the 



spatial metric, Uij := dtgij by 



Ma = — tin. 
13 2a 3 



(35) 



Because My is trace free [Eqs. (13) and (15)-(16)], we 
require Uij to be trace free. 

The conformal thin sandwich equations allow control 
of certain time-derivatives in the subsequent evolution 
of the constructed initial data. If the lapse a and shift 
(3 l from the initial data are used in the evolution, for 
instance, then the trace-free part of d t will be propor- 
tional to Uij. Therefore (see Refs. [27, 30]) 







(36a) 



is a preferred choice for initial data sets that begin nearly 
in equilibrium, such as binary black holes quasi-circular 
orbits. 

The evolution equation for K can be used to derive an 
elliptic equation for the conformal lapse a (or, cquiva- 
lcntly, for aijj). Upon specification of 



d t K = 0, 



(36b) 



this fifth elliptic equation is to be solved for a simulta- 
neously with Eqs. (17) and (18), cf. [27, 30]. 



G 



Our numerical code uses the conformal factor ip, the 
shift /3 1 , and the product of lapse and conformal factor 
aip = aip 7 as independent variables, in order to simplify 
the equation for dtK. Thus, the actual equations being 
solved take the form 



V 2 f 



-Rip 



1 

12 



(37a) 



(Lpy 



o = 



2(aip) 

i\P 
.2(aip) 

V 2 (aip) - (aip) 

iP 5 (d t K~/3 k d k K) 



(37b) 



with 



i\P 
2mp 



(37c) 



(37d) 



These equations can be solved only after 

1. specifying the remaining free data: i.e., the confor- 
mal metric g^ and the trace of the extrinsic curva- 
ture K (we chose already Uij = and d t K = 0), 

2. choosing an inner boundary S which excises the 
black holes' singularities, and also an outer bound- 
ary B, and 

3. choosing boundary conditions for ip, aip, and (3 l on 
B and S. 

The initial data is required to be asymptotically flat, 
and the outer boundary B is placed at infinity 7 . If 
is asymptotically flat, the outer boundary conditions are 
then 





aip 



1 on B, 
1 on B, 
(ft x r) 



+ dor 1 on B. 



(38a) 
(38b) 
(38c) 



Here r l is the coordinate position vector. The shift 
boundary condition consists of a rotation (parametrized 
by the orbital angular velocity fJo) and an expansion 
(parametrized by do); the initial radial velocity is neces- 
sary for reducing orbital eccentricity in binary-black-holc 
initial data [49]. 

The inner boundary condition on the conformal fac- 
tor ip ensures that the excision surfaces S are apparent 
horizons [27]: 



8a 



0" 



(L/3) i:; - - 



i°3 



-KiP 3 on S. 
6 



(39) 



Here s l := ip 2 s l , s 1 is unit vector normal to S, and hi 
gij — §iSj is the induced conformal 2-metric on S. 
The inner boundary condition on the shift is 



P = as 1 - Q r C on S, 



(40) 



where £ l Si = 0. The first term on the right-hand-side 
ensures that the apparent horizons are initially at rest; 
the tangential term determines the black hole's spin [27- 
29]. 

References [27-29] chose the sign of the last term in 
Eq. (40) such that positive values of f2 r counteract the 
spin of the corotating holes that are obtained with Q r = 
0. Here, we are interested in large spins, and we reverse 
the sign of the last term in Eq. (40) so that positive, 
increasing f2 r results in increasing spins. 

Two sets of choices for g^, K, S, and the boundary 
condition for aip on S are discussed in the next sub- 
sections. Each set of choices will be used to construct 
binary-black-hole initial data in Sec. IV. 



1. Conformal flatness & maximal slicing (CFMS) 

The simplest choice for g^ is a flat metric, 

9ij = fij- (41) 

This choice has been used almost exclusively in the pre- 
vious formulations of binary-black-hole initial data. 

The simplest choice for K, also commonly used in prior 
formulations of binary-black-holc initial data, is maximal 
slicing, i.e. 



K = 0. 



(42) 



Also for simplicity, we choose to make the excision sur- 
face S consist of coordinate spheres: 



S = (J S a , 



(43) 



where iS a are surfaces of constant Euclidean distance r exc 
about the center of each excised hole, and n — 1 or 2 is 
the number of black holes present in the initial data. 

The boundary condition for the lapse on S determines 
the temporal gauge; we adopt the condition given in 
Eq. (59a) of Ref. [28]: 



d 
dr a 



(aip) = on S a , 



(44) 



where r a is the Euclidean distance from the center of hole 
a. This type of initial data is used in Refs. [49-51]. 



2. Superposed Kerr Schild (SKS) 



7 In practice, B is a sphere with radius > 10 9 times the coordinate 
radius of the black-hole horizons. 



Single black holes with angular [52, 53] or linear [54] 
momentum do not admit conformally-flat spatial slicings; 
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therefore, conformal flatness [Eq. (41)] is necessarily defi- 
cient. This has motivated investigations of binary-black- 
hole initial data whose free data have stronger physical 
motivation, e.g. Refs. [37, 38, 55-61]. 

In this subsection, we consider conformally-curved 
data that are in the same spirit as the SKS data of 
Refs. [37, 38] although here i) we apply the idea to the 
QE-XCTS formalism, and ii) as discussed below, our free 
data is very nearly conformally-flat and maximally-sliced 
everywhere except in the vicinity of the black holes. 

The choices we make here generalize the conformally- 
curved data in chapter 6 of Ref. [39] to nonzero spins. 
Specifically, the free data and lapse boundary condition 
will be chosen so that the conformal geometry near each 
hole's horizon is that of a boosted, spinning, Kerr-Schild 
black hole. The conformal metric gij and the mean cur- 
vature K take the form 

n 

hi ■= /« + X>-^(fl5-/y), (45) 
n 

K := Y. e ' rl ' wlK - ( 46 ) 



Here gfj and K a are the spatial metric and mean cur- 
vature, respectively, of a boosted, spinning Kerr-Schild 
black hole with mass M a , spin S a , and speed v a . 

Far from each hole's horizon, the conformal metric is 
very nearly flat; this prevents the conformal factor from 
diverging on the outer boundary [39]. The parameter 
w a is a weighting factor that determines how quickly the 
curved parts of the conformal data decay with Euclidean 
distance r a (a = 1,2, ...) from hole a; in this paper, the 
weight factor w a is chosen to be larger than the size scale 
of hole a but smaller than the distance d to the com- 
panion hole (if any): M a < w a ;$ d a . This is simi- 
lar to the "attenuated" superposed-Kerr-Schild data of 
Refs. [38, 62], except that here the weighting functions 
are Gaussians which vanish far from the holes, while in 
Refs. [38, 62] the weighting functions go to unity far from 
the holes. 

The excision surfaces S a are not coordinate spheres 
unless S a — and v a — 0. Instead they are deformed in 
two ways, i) They are distorted so that they are surfaces 
of constant Kerr radius rjcerrj i-c 



x 2 + y 2 



1 Kerr 



S a /Ma 



1 



(47) 



where x, y, and z are Cartesian coordinates on the S. 
Then, ii) the excision surfaces are Lorentz-contracted 
along the direction of the boost. 

The boundary condition for the lapse a on S a is a 
Dirichlct condition that causes a (and, consequently, the 
temporal gauge) in the vicinity of each hole to be nearly 
that of the corresponding Kerr-Schild spacetime, i.e. 



1 



-r a /w' c c 



(a a - 1) on S a , 



(48) 



where a a is the lapse corresponding to the Kerr-Schild 
spacetime a. 



III. SINGLE-BLACK-HOLE INITIAL DATA 
WITH NEARLY-EXTREMAL SPINS 

In this section, we examine to which extent the for- 
malisms presented in Sec. II can generate single black 
hole initial data with nearly-extremal spin. We consider 
first Bowen-York initial data and then conformally-flat 
quasi-equilibrium data. Since superposed-Kerr-Schild 
data can represent single Kerr black holes exactly, there is 
no need to investigate single-hole superposed-Kerr-Schild 
data. In Sec. IV, we will both consider conformally-flat 
and superposed-Kerr-Schild data for binary black holes. 

To orient the reader, the initial data sets constructed 
in this section, as well as the binary-black-hole data sets 
constructed in Sec. IV, are summarized in Table I. 

Unless noted otherwise, all spins presented in this sec- 
tion are measured using the approximate-Killing-vector 
spin xakv described in Appendix A. Therefore, the sub- 
script "AKV" in xakv will be suppressed for simplicity. 



A. Bowen-York (puncture) initial data 

As discussed in Sec. II B, for a single spinning black 
hole at rest, puncture initial data is identical to inversion- 
symmetric initial data. Such solutions have been exam- 
ined in the past (e.g. [48, 63]), and additional results were 
obtained (partly in parallel to this work) in the study by 
Daiii, Lousto, and Zlochower al [25]. 

We revisit this topic here to determine the maximum 
possible spin of Bowen-York (BY) initial data more accu- 
rately than before, to establish the power-law coefficients 
for the approach to these limits with increasing spin pa- 
rameter S 1 , and to present new results about the geomet- 
ric structure of Bowen-York initial data with very large 
spin parameter. 

We solve Eq. (30) with the pseudo-spectral elliptic 
solver described in Ref. [64]. The singular point of u 
at the origin is covered by a small rectangular block ex- 
tending from ±10 _4 TOp along each coordinate axis. This 
block overlaps four concentric spherical shells with radii 
of the boundaries at 8 • 10~ 5 m p , 0.005m p , 0.3m p , 50m p , 
and 10 9 m p . The equations are solved at several different 
resolutions, with the highest resolution using 20 3 basis- 
functions in the cube, L = 18 in the spheres and 26 
and 19 radial basis-functions in the inner and outer two 
spherical shells, respectively 

Because of the axisymmetry of the data-set, the ro- 
tational Killing vector of the apparent horizon is sim- 
ply £?0. The integral for the quasilocal spin, Eq. (Al) 
turns out to be independent of ip and can be evaluated 
analytically with a result equal to the spin-parameter, 
S. Thus we can use this initial data set to check how 
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Label 


Section 


Figures 


n 


d 


fio 


do x 10 4 


fir or S/rrip 


S 


|XAKV| 


M iTT 


M 


-Eadm 


BY-Single 


III A 


1-5, 8, 19 


1 


- 


- 


- 


0.01 < S/ml < 10 4 


- 










CFMS-Single 


IIIB 


6-8, 19 


1 


- 


- 


- 


< fir < 0.191 


- 










CFMS 


IV A 


9, 13 


2 


32 


0.007985 





< fi r < 0.1615 


- 










SKS-0.0 


IV B 


11, 13 


2 


32 


0.006787 





< fi r < 0.24 













SKS-0.5 


IV B 


11, 13 


2 


32 


0.006787 





< fi r < 0.27 


0.5 










SKS-0.93 


IV B 


11-13 


2 


32 


0.006787 





< fir < 0.35 


0.93 










SKS-0.99 


IV B 


10-13 


2 


32 


0.007002 


3.332 


0.28 < fi r < 0.39 


0.99 










SKS-0.93-E0 


VB 


14 


2 


32 


0.006787 





0.28 


0.93 


0.9278 


0.9371 


1.131 


2.243 


SKS-0.93-E1 


VB 


14 


2 


32 


0.007 





0.28 


0.93 


0.9284 


0.9375 


1.132 


2.247 


SKS-0.93-E2 


VB 


14 


2 


32 


0.006977 


3.084 


0.28 


0.93 


0.9275 


0.9395 


1.134 


2.249 


SKS-0.93-E3 


VC 


10-11, 13-16, 19 


2 


32 


0.007002 


3.332 


0.28 


0.93 


0.9275 


0.9397 


1.134 


2.250 


SKS-HeadOn 


VD 


10-11, 13, 17-19 


2 


100 








0.3418 


0.97 


0.9701 


0.8943 


1.135 


2.257 



TABLE I: Summary of the initial data sets constructed in this paper. The first row (BY-Single) represents Bowen-York initial 
data for single black holes of various spins. The next two rows (CFMS-Single and CFMS) are quasi-equilibrium, conformally- 
flat, maximally-sliced initial data for single and binary spinning black holes, respectively. All other data sets employ superposed 
Kerr-Schild quasi-equilibrium data with the second block of rows representing families of initial data sets for various spins and 
the last block of rows representing individual data sets to be evolved. The data sets SKS-0.93-E0 to SKS-0.93-E3 demonstrate 
eccentricity removal, and SKS-HeadOn is used in a head-on evolution. The first block of columns gives the label used for each 
data set, and the relevant section of this paper devoted to it. The next block of columns lists the most important parameters 
entering the initial data. The last block of columns lists some properties of those data sets that we evolve in Sec. V. 
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FIG. 1: Convergence test for a single puncture black hole 
with a very large spin parameter S/rrip = 10000. Plotted 
are results vs. resolution N, which is the total number of 
basis-functions. The solid lines show the relative differences 
of three angular momentum measures to the analytically ex- 
pected value 10000. The dashed lines show differences from 
the next-higher resolution of two dimensionless quantities for 
which no analytic answer is available. 



well our spin-diagnostics and our ADM angular momen- 
tum diagnostic works (recall that Jadm is also equal to 
the spin-parameter S). This comparison is performed 



in Fig. 1, which shows relative differences between the 
numerically extracted values for the approximatc-Killing- 
vector (AKV) spin, the coordinate spin (defined with the 
AKV spin in Appendix A), and the ADM angular mo- 
mentum Jadm relative to the expected answer, S. The 
figure also shows differences between neighboring resolu- 
tions for the two quantities of interest below, S/M 2 = \ 
and Syi? ADM = Jadm/^adm = £ J- 



Figure 1 seems to show exponential convergence with 
increased resolution N. Since puncture data is only C 2 at 
the puncture, one would rather expect polynomial con- 
vergence. The effect of the non-smoothness at the punc- 
ture is mitigated by choosing a very high resolution close 



to the puncture (a small cube with sides ±10 



with 



20 3 basis- functions). Therefore, for the resolutions con- 
sidered in Fig. 1, the truncation error is dominated by 
the solution away from the puncture, and exponential 
convergence is visible. If we used infinite-precision arith- 
metic and were pushing toward higher resolution than 
shown in Fig. 1, then we would expect to eventually sec 
polynomial convergence dominated by the cube covering 
the puncture. 



Next, we construct a series of initial data sets with 
increasing spin-parameter S, and compute x, £j, and £ 
for each initial data set. The results are plotted in Fig. 2 
and confirm earlier results [26, 63]. In addition, the inset 
shows that the asymptotic values Xmax = 0.9837 and 
£u max = 0.928200 arc approached as power-laws in the 
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10 100 1000 10000 



S/nf 




0.001 



2r/m = r/Rinv 

p 



FIG. 2: Properties of single, spinning puncture black holes 
with spin-parameter S and puncture mass m p . The di- 
mensionless spin \ := S/M 2 , ADM angular momentum 
£.j '■= Jadm/-E1dmi an d spin-extremality parameter C, '■= 
S/ (2M^ r ) are plotted against the spin parameter S/m 2 ,. The 
horizon mass M is related to the spin S and irreducible mass 
M irl in Eq. (2). 



spin-parameter, 



Xmax ~ X OC 



max £j ®£ I o 



5' 



-0.75 



-1.4 



(49) 
(50) 



The exponents of these power-laws are computed here for 
the first time. 

To confirm that the apparent horizon is indeed at 
r = i?; nv , we ran our apparent horizon finder on the 
high-spin puncture initial data sets. The horizon finder 
had great difficulty converging, and the reason for this 
becomes clear from Fig. 3. The main panel of this fig- 
ure shows the area of spheres with coordinate radius r. 
The area is minimal at r = m p /2, as it must be, since 
rrip/2 = i?i nv is the radius of the inversion sphere. How- 
ever, the area is almost constant over a wide range in 
r — for S/rrip = 10000 over about two decades in either 
direction: 0.01 < r/R im < 100. Thus, the Einstein- 
Rosen bridge (the throat) connecting the two asymptot- 
ically flat universes lengthens as the spin increases, giv- 
ing rise to an ever-lengthening cylinder. If this were a 
perfect cylinder, then the expansion would be zero for 
any r = const cross-section. Because the geometry is 
not perfectly cylindrical, the expansion vanishes only for 
r = m v /2 = i?i nv , but remains very small even a sig- 
nificant distance away from r = m p /2 = i?i nv . This is 



FIG. 3: Properties of coordinate spheres with radius r for 
high-spin puncture initial data. Main panel: Area of these 
spheres. Inset: residual of the apparent horizon equation on 
these spheres. The area is almost constant over several orders 
of magnitude in r. The apparent-horizon-residual vanishes at 
r = Rinv, but is very small over a wide range of r. 



shown in the inset, which plots the residual of the appar- 
ent horizon finder at different radii. 

With the lengthening of the throat, the interval in r 
with small expansion lengthens, and the value of the ex- 
pansion within this interval reduces. Both effects make 
it harder for the apparent horizon finder to converge. In 
Fig. 2, we have used our knowledge of the location of 
the apparent horizon to set tah = rnp/2, rather than 
to find this surface numerically. Without this knowl- 
edge, which arises due to the identification of puncture 
data and inversion symmetric data, computation of Fig. 2 
would have been significantly harder, perhaps impossible. 

Let us assume for the moment that the solution ip( r ) = 
-j^ + l + u(r) is spherically symmetric (we give numerical 
evidence below that this is indeed a good approximation). 
Because gij = ip 4 fij, the area of coordinate spheres is 
then given by 



A(r) = ATiij} 2 {r)r. 



(51) 



In the throat region, where A(r) w const, the conformal 
factor must therefore behave like l/^/r 1 as also argued 
independently by Dain, Lousto, and Zlochower [25]. 

To extend on Dain et aL's analysis, let us substitute 
Eq. (23) into Eq. (24) to obtain the well-known equation 



9S 2 sin 2 9 
4r 6 



(52) 



where 9 is the angle between the spin-direction and the 
point x l . Continuing to assume that ip is approximately 
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spherically symmetric, we can replace the factor sin 2 by 
its angular average (47r) _1 J sin 2 6 dVl = 2/3, and obtain 



"I I I 



d 2 ip + 2dip 
dr 2 r dr 



3S 2 
2r 



(53) 



Here, we introduced an overbar ip to distinguish the 
spherically symmetric solution ip(r) of Eq. (53) from the 
full solution ip(x l ) of puncture/inversion-symmetric ini- 
tial data. Following Dain et al. [25] we assume that the 
conformal factor behaves as a power-law (ip(r) = Ar a ) 
and substitute this into Eq. (53). We find that Eq. (53) 
determines the power-law exponent a = —1/2 and the 
overall amplitude A = (6S 2 ) 1 / 8 , so that 



ip{r) 



(6S 2 



96 1 / 8 f ^ 

. nip 



1/4 



-1/2 



(54) 



In Eq. (54), we chose the scaling S/m 2 which is com- 
monly used in the puncture-data literature, but kept 
r/Rinv to emphasize the inversion symmetry of the data 
in our figures (in a log-plot using r/i? inv , the solution will 
appear symmetric, see e.g. Fig. 3). While ip(r) solves 
the spherically symmetric Eq. (53) exactly, it must de- 
viate from ip(x l ) for sufficiently large r because ip — > 
as r — > oo, whereas ip — > 1. The deviation will become 

significant when ip ~ 1, i.e. at radius r x ~ ^Js/m 2 . Be- 
cause of inversion symmetry, this implies a lower bound 
of validity at l/r x , so that Eq. (54) holds for 



S 



-1/2 



< 



< 



S 



1/2 



The circumference of the cylindrical throat is 



C = 2ixip{rfr = 2tt96 1/4 




(55) 



(56) 



and its length is 



(S/m 



{S/mlY 



2,1/2 



^ 2 {r) dr = 96 1 / 4 ,/A m ( 4r 
rrip \ m 



Therefore, the ratio of length to circumference, 



C 2tt 



(57) 



(58) 



grows without bound as S/m 2 becomes large, albeit very 
slowly. The scaling with [S/m 2 ) 1 ! 2 in Eqs. (55)-(57) 
might seem somewhat surprising. However, in the large 
spin limit, S/M 2 is just a constant close to unity (namely 
Xmax = 0.9837). Therefore, S 1 / 2 w M, i.e. the scaling 
S 1/>2 is effectively merely a scaling with mass. 

Figure 4 shows the conformal factor ip, the "puncture 
function" u, and the estimate ip of Eq. (54) for three 
different values of S/m 2 . There are several noteworthy 
features in this figure. First, both ip and u show clearly 
three different regimes: 



1000 



¥ 

— u 




S/m =10000 

p 

1000 
100 



0.001 0.01 



III I I 

0.1 1 10 
2r/m = r/R. 

p inv 



FIG. 4: Solutions of high-spin puncture initial data. Plotted 
are the conformal factor ip and puncture function u in the 
equatorial plane as a function of radius r. Furthermore, the 
approximate solution ip is included, with solid circles denot- 
ing the range of validity of this approximation, cf. Eq. (55). 
Three curves each are plotted, corresponding from top to bot- 
tom to S/m 2 v = 10000, 1000, 100. 



• For large r, ip ss 1 and u a 1/r, 
asymptotically-flat end. 



This is the upper 



• For intermediate r, ip oc 1/v^* an d u l/v 7 ^ 7 - 
This is the cylindrical geometry extending sym- 
metrically around the throat. This region becomes 
more pronounced as S increases. 



• For small r, ip oc 1/r and u ? 
lower asymptotically-flat end. 



const. This is the 



Figure 4 also plots the approximate solution ip [cf. 
Eq. (54)] for its range of validity_ [given by Eq. (55)]. 
Note that slope and amplitude of ip fit very well the nu- 
merical solution ip. In fact, the agreement is much better 
than with u. 

One could also have started the calculation that led to 
Eq. (54) with Eq. (30). Assuming spherical symmetry 
and assuming that u ^> m p /{2r) + 1, we would have de- 
rived Eq. (53), but with ip replaced by u. We would then 
have found the approximate behavior Eq. (54) for u. The 
disadvantage of this approach is the need for additional 
approximations, which reduce the accuracy of the result. 
From Fig. 4 we see that, in the throat region, the dotted 
lines representing ip are close to the dashed lines of u. 
But the agreement between ip and ip is certainly better. 

Finally, we note that the limits of validity of ip 
[Eq. (55)] match very nicely the points where the nu- 
merical ip diverges from ip. 
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p inv 



FIG. 5: Angular decomposition of the conformal factor 
if)(r,9,<f>) for single black hole puncture data. 



To close this section, we present numerical evidence 
that indeed ip is approximately spherically symmet- 
ric, the assumption that entered into our derivation of 
Eq. (54). We decompose the conformal factor of the nu- 
merical puncture data solutions into spherical harmonics, 

oo I 

iP(r,9,<t>)=J2 E feW^W). ( 59 ) 

1=0 m=-l 

and plot in Fig. 5 the sizes of the I ^ modes rela- 
tive to the spherically symmetric mode V'oo- Because of 
the symmetries of the problem, the only non-zero modes 
have m = and even In the throat region, the largest 
non-sphcrically symmetric mode ?p20 is about a factor of 
65 smaller than the spherically symmetric mode. With 
increasing /, ipi m decays very rapidly. Also, in both 
asymptotically flat ends, the non-spherically symmetric 
modes decay more rapidly than the I = mode, as ex- 
pected for asymptotically flat data. This figure again 
shows nicely the inversion symmetry of the data, under 
f 7 R'mv (r 7 'i?inv) _1 • Given the simple structure of the 
higher modes, it should be possible to extend the ana- 
lytical analysis of the throat to include the non-spherical 
contributions. To do so, one would expand ip as a se- 
ries in Legcndrc polynomials in 6; the ip~ 7 -term on the 
right hand side of Eq. (52) would result in a set of ordi- 
nary differential equations for those coefficients. In the 
throat region, the radial behavior of each mode should be 
oc 1/y/r, and the ordinary differential equations should 
simplify to algebraic relations. 




0.05 0.1 0.15 0.2 

Q. 

r 

FIG. 6: Conformally-fiat, maximally-sliced, quasiequilibrium 
initial data sets with a single, spinning black hole. We plot the 
horizon mass M, irreducible mass Mi rr , and the (approximate- 
Killing-vector) spin S against the rotation parameter Q r [cf. 
Eq. (40)]. Only Q r is varied in this figure; all other parameters 
are held fixed. The upper and lower points with the same fir- 
are obtained numerically by choosing different initial guesses. 
The inset shows a close-up view of the turning point, which 
occurs at Sl r ~ 0.191. 



B. Quasi-equilibrium 
extended-conformal-thin-sandwich data 

We have seen in Sec. Ill A that puncture initial data 
for single, spinning black holes can be constructed for 
holes with initial spins of \ < 0.9837. In this section, 
we address the analogous question for excision black- 
hole initial data: how rapid can the initial spin be for a 
single, spinning black hole constructed using quasiequi- 
librium, cxtendcd-conformal-thin-sandwich (QE-XCTS) 
initial data? 

As noted previously, if the free data gij and K are 
chosen to agree with the analytic values for a Kerr black 
hole, g£ crr and K KcTI , then the QE-XCTS initial data can 
exactly represent a single Kerr black hole. In this case, 
X = 1 is obtained trivially by choosing S = M 2 = 1, 
where M and S are the mass and spin, respectively, of 
the Kerr black hole described by the conformal metric. 

Setting aside this trivial solution, we construct 
conformally-flat, maximally-sliced (CFMS) data for a 
single, spinning hole. We construct a family of QE- 
XCTS initial data sets for single spinning black holes 
by numerically solving the XCTS equations [in the form 
stated in Eqs. (37a)-(37c)] using the same spectral ellip- 
tic solver [64] as in Sec. Ill A. The free data are given by 
Eqs. (41)-(42) and by Eqs. (36a)-(36b). 

On the outer boundary B, we impose Eqs. (38a)-(38c). 
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So that the coordinates are asymptotically inertial, we 
choose Qo = a o = in Eq. (38c). 

We excise a coordinate sphere of radius r cxc about the 
origin, where 



0.85949977 



(60) 



is chosen such that for zero spin M = 1. On this inner 
boundary S, we impose Eqs. (39)-(40) and Eq. (44). The 
spin is determined by Eq. (40): first, the vector £ 8 is cho- 
sen to be the coordinate rotation vector 3$, making the 
spin point along the positive z axis; then, the rotation pa- 
rameter fl r is varied while the other parameters are held 
fixed. The spin is measured on the apparent horizon us- 
ing the approximate-Killing- vector spin (Appendix A); 
because in this case the space is axisymmetric, the "ap- 
proximate" Killing vector reduces to the corresponding 
exact rotational Killing vector. 

Figure 6 show how the mass M and AKV spin S 
depend on £l r . At £l r = 0, we find the spherically- 
symmetric solution with S = and M; rr = M = 1 (the 
mass is proportional to the excision radius, and Eq. (60) 
sets it to unity). Using this spherically-symmetric so- 
lution as an initial guess for the elliptic solver, we find 
solutions for increasing Q r with spin increasing initially 
linearly with Q r and with approximately constant mass. 
Beyond some critical n rjCr it, the elliptic solver fails to 
converge, and close to this point, all quantities vary in 
proportion to yfV^rit — These symptoms indicate 
a critical point where the solutions "turn over" and con- 
tinue towards smaller f2 r . Analogous non- unique solu- 
tions of the XCTS equations have been discovered be- 
fore in Ref. [34]. To construct solutions along the up- 
per branch, one must choose a sufficiently close initial 
guess for the elliptic solver; we follow the steps outlined 
in Ref. [34] and are able to find solutions along the upper 
branch for a wide range of f2 r < f2 rjC rit- As Fig. 6 shows, 
mass and spin of the horizon in solutions along the up- 
per branch increase with decreasing f2 r , analogous to the 
findings in [34, 35]. 

Figure 7 shows the dependence of % = S/M 2 , ej = 
4dm/-Badmi an d C = S/ (2M£ r ) on O r . The curves 
reflect again the non-unique solutions. The dimension- 
less spin \ increases continuously along the lower branch, 
and reaches x ~ 0-85 at the critical point. As fl r is de- 
creased along the upper branch, x continues to increase, 
eventually reaching values larger than 0.99. It appears 
X continues to increase as il r — > 0. To find the limiting 
value, consider that the behavior of the extremality pa- 
rameter £ in the inset of Fig. 7. Assuming that £ can 
be extrapolated to O r — > 0, we find a limiting value of 
C ~ 0.88. By Eq. (9), this implies a maximal value of 
X ~ 0.992. 

In Figs. 6-7, the data sets on the lower branch ap- 
pear to be physically reasonable. For spins x ^5 0.85, 
the mass M is nearly constant, and the dimension- 
less spin x increases linearly with fl r . Furthermore, as 
f2 r — > the lower branch continuously approaches the ex- 
act Schwarzschild spacetime (see [28] ) . The upper branch 
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FIG. 7: Conformally-flat, maximally-sliced quasiequilibrium 
initial data sets with a single spinning black hole: The di- 
mensionless spin \, dimensionless ADM angular momentum 
ej, and spin-extremality parameter £ plotted against fi r [cf. 
Eq. (40)]. Only Q r is varied in this figure; all other parame- 
ters are held fixed. The inset enlarges the area in the upper 
left corner; we are able to generate data sets with \ > 0.99, 
whereas the largest spin obtainable on the lower branch is 
X ~ 0.85. 



appears to be physically less reasonable; for instance, 
the spin x increases for decreasing horizon frequency f2 r . 
Comparing Figs. 2 and 7, we see that the QE-XCTS data 
leads to somewhat larger values of x an d £ J relative to 
puncture data. However, the values are not too different, 
and similar trends remain. For instance, x is much closer 
to unity than ej. 

To investigate differences or similarities between punc- 
ture data and QE-XCTS data further, we compute em- 
bedding diagrams of the equatorial planes of these data 
sets. The initial data for single black holes have rota- 
tional symmetry about the z-axis, so the metric (12) on 
the initial data hypcrsurface, when restricted to the equa- 
torial plane, can be written as 



ds 2 = V/ {dr 2 + r 2 d(j> 2 ) 



(61) 



where r and <j) are the usual polar coordinates. This 
metric is now required to equal the induced metric on 
the 2-D surface given by Z = Z(R) embedded in a 3-D 
Euclidean space with line-clement 



°^Euclidcan ~ dR 



R" 



(62) 



Setting dZ = ^|dR, we obtain the induced metric on the 
Z = Z(R) surface 



ds 2 = 



dZ 
dR 



dR 2 + R 2 d6 2 . 



(63) 
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FIG. 8: Embedding diagrams for puncture and quasiequilib- 
rium initial data. Plotted is the embedding height Z as a func- 
tion of the embedding radius R, both scaled by the mass M. 
For quasiequilibrium data (dashed lines), Z=0 at r = r cxc ; for 
puncture data (solid lines), Z=0 at r = Ri nv . The thin solid 
purple curve represents the embedding of a plane through a 
Schwarzschild black hole in Schwarzschild slicing. 



10 while m p = 1 is held constant, we find from Eq. (57) 
that C/S 1 / 2 should increase by 

Ofil/4 

AC/S 1/2 = — — In 10 « 3.60, (67) 

where the factor 1/2 arises because i?i nv = m p /2 = 0.5. 
The embedding diagram shows only the top half of the 
throat, and 5 1 / 2 m M [cf. the discussion after Eq. (58)]. 
Therefore in Fig. 8 the S = 100, 1000, 10000 lines for BY 
(puncture) data should be spaced by AZ/M w 1.80 for 
large R/M . This indeed is the case. 

The CFMS datasets appear to scale proportionally to 
y/S, which is similar to the puncture data's behavior. 
Furthermore, the CFMS initial data sets also develop a 
lengthening throat as S becomes large (the effect is not 
as pronounced as for puncture data, owing to the smaller 
maximal S we achieved.) Thus it appears that large spin 
CFMS data might be similar to large spin puncture data. 
However, the throats of the QE-XCTS data show a bulge 
near the bottom, because for these data sets R actually 
decreases with r in the immediate vicinity of r cxc . This 
is unlike the puncture data, which very clearly exhibit 
cylindrical throats, consistent with the discussion leading 
to (58). 



Equating Eqs. (61) and (63), we find 
R = ip 2 r 

and 



1 



dZ 
dR 



dR 2 = ip A dr 2 . 



(64) 



(65) 



Combining (65) and (64) results in 



dZ 
dr 



-Artj} 2 —— ( tp + r-f- 
dr \ dr 



(66) 



Since the pseudo-spectral elliptic solver gives ip as a func- 
tion of r, Eqs. (64) and (66) allow us to solve for the em- 
bedding radius R and the embedding height Z in terms 
of r. 

Figure 8 shows embedding diagrams for three sets of 
QE-XCTS and puncture data. We have set Z=0 at 
T = r cxc for QE-XCTS data and at r = i?mv for punc- 
ture data. This figure also contains the embedding of 
a plane through Schwarzschild in Schwarzschild coordi- 
nates (i.e. the 5 = limit of BY puncture data), given by 
R/M = Z 2 /{8M 2 ) + 2. Both puncture data and CFMS 
data exhibit a lengthening throat with increasing spin 
S/M 2 . For puncture data, this lengthening can be de- 
duced from the analytical results in Sec. Ill A: as the spin 
parameter S of the puncture data increases by a factor of 



IV. BINARY-BLACK-HOLE INITIAL DATA 
WITH NEARLY-EXTREMAL SPINS 

In this section, we construct binary-black-hole initial 
data with rapid spins, confining our attention to the spe- 
cial case of spins aligned with the orbital angular mo- 
mentum. In the limit of large separation, binary-black- 
hole puncture initial data will behave like two individual 
puncture initial data sets. Specifically, we expect that 
it should be possible to construct puncture binary-black- 
holc initial data with initial spins x(t — 0) ^ 0.98, but 
the spins will rapidly drop to x < 0.93 as the black holes 
settle down. For this reason, and also because puncture 
data is not well-suited to our pseudospectral evolution 
code, we will restrict our attention to binary black holes 
constructed with the QE-XCTS approach. 

As laid out in Table I, we first construct a family (la- 
belled CFMS) of standard conformally-flat initial data 
on maximal slices; then, we turn our attention to families 
(labelled SKS) of superposed Kerr-Schild initial data. Fi- 
nally, we construct a few individual SKS initial-data sets 
which we evolve in Sec. V. All of the data sets represent 
equal-mass, equal-spin black holes with spins parallel to 
the orbital angular momentum. 

In this section, unless otherwise indicated, all dimen- 
sionlcss spins are the approximatc-Killing-vector spin 
Xakv (Appendix A), and the subscript "AKV" will be 
suppressed for simplicity. 
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FIG. 9: Main panel: Dimensionless spin x [Eq. (1)] and spin- 
extremality parameter ( [Eq. (8)] for the family CFMS of 
spinning binary-black-hole initial data. Inset: Enlargement 
of x toward the end of the upper branch, with circles de- 
noting the individual initial data sets that were constructed. 
Compare with Fig. 7. 



A. Conformally flat, maximal slicing data (CFMS) 

To construct conformally-flat binary-black-hole data, 
we solve the same equations and boundary conditions as 
for the single-black-hole described in Sec. IIIB, 

with the main difference being that we excise two spheres 
with radius r cxc [cf. Eq. (60)] with centers on the x- 
axis at x = ±d/2. The initial spins of the holes are 
set by adjusting Q r , just as in the single- hole case. The 
parameters f2o and clq in the outer boundary condition 
on the shift [Eq. 38c] determine the initial angular and 
radial motion of the holes, which in turn determine the 
initial eccentricity e of the orbit. We set flo = ilo e z, 
where e z is a unit vector that points along the positive z 
axis. For the CFMS family of data sets considered here, 
we use values for fio and do that should result in closed, 
fairly circular orbits, since our choices of £lo and do lead 
to data sets that approximately satisfy the Komar-mass 
condition -Eadm = Mr (cf. [29]). Specifically, on the 
lower branch of the resulting non-unique family of initial 
data, 



\E AB m-M k \ 



< 1%, 



M 



(68) 



where the Komar mass is defined by (e.g., Eq. (35) of 
Ref. [29]) 



dA. 



(69) 



(On the upper branch, £adm and Mk differ by up to 
3%.) 

As the rotation parameter f2 r is varied (with the coor- 
dinate separation d held fixed), we find that the CFMS- 



family of binary-black-holc initial data behaves qualita- 
tively similarly to the analogous single-black-hole initial 
data discussed in Sec. Ill B. There is a maximal fi riC rit 
such that no solutions can be found for fi r > f2 r)Cr it; for 
values of O r below f^crit, two solutions exist. Figure 9 
plots the dimensionless spin \ an d the spin-extremality 
parameter £ against £l r for this family of initial data. We 
only show values for one of the holes, since the masses 
and spins are equal. Spins larger than \ ss 0.85 appear 
on the upper branch. The highest spin we have been able 
to construct is larger than x = 0.97. 



B. Superposed-Kerr-Schild data 

In this section, we solve the same equations and bound- 
ary conditions as in the conformally flat case, except that 
we use SKS free data (Sec. II C 2) instead of conformally- 
flat free data. To construct the individual Kerr-Schild 
data, we need to choose for each black hole the coordi- 
nate location of its center, its conformal mass M, confor- 
mal spin S, and its boost-velocity. We center the black 
holes on the x-axis at x = ±d/2, use the same mass 
M = 1 for both black holes, and set the boost velocity 
to (0, ±d fio/2, 0). The conformal spins are always equal 
and arc aligned with the orbital angular momentum of 
the holes. 

In contrast to the CFMS data, there are now two pa- 
rameters that influence the black holes' spins: i) the 
rotation parameter Cl r in Eq. (40), and ii) the confor- 
mal spin S. For concreteness, we choose to construct 
data for four different values of the conformal spin: 
S/M 2 = 0,0.5,0.93, and 0.99. For each choice, we con- 
struct a family of initial data sets for different values of 
Q r , which we label as SKS-0.0, SKS-0.5, SKS-0.93, and 
SKS-0.99 respectively. 

Other choices that went into the construction of the 
SKS initial data sets are as follows: 

• The excision boundaries are chosen to be the co- 
ordinate locations of the horizons of the individual 
Kerr-Schild metrics, i.e. they are surfaces of con- 
stant Kcrr-radius 



d "' = r+ := M + VAP - S 2 , 



(70) 



length-contracted by the Lorentz-factor appropri- 
ate for the boost velocity of each black hole. This 
length-contraction accounts for the tangential mo- 
tion of the hole but neglects the much smaller radial 
motion. 



• When superposing the individual Kerr-Schild met- 
rics, we use a damping length scale w — 10r^,° c lr [cf. 
Eqs. (45) and (46)], except for the SKS-0.99 family, 
which uses w = d/3. 

• The orbital frequency Qq and radial expansion do 
are held fixed along each family. We expect that 
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10: Convergence of the spectral elliptic solver. Left 
The residual constraint violation as a function of the 
total number of grid-points N when running the elliptic solver 
at several different resolutions. Right panel: Convergence of 
the black hole dimensionless spin \ [Eq. (1)] with increasing 
resolution Lah of the apparent horizon finder, applied to the 
highest-resolution initial data set of the left panel. The three 
curves in each panel represent three different initial data sets: 
One from the family SKS-0.99, as well as the two initial data 
sets that are evolved in Sec. V. 



our choices for f^o and do will lead to bounded, 
fairly circular orbits, since 



I-Badm — Mk\ 
-Eadm 



< 3%. 



(71) 



In Sec. VB we reduce the orbital eccentricity for 
one data set in the family SKS-0.93. 

We again solve the XCTS equations using the spectral 
elliptic solver of Ref. [64] ; the families of SKS initial data 
sets that we construct are summarized in Table I. The 
elliptic solver needs some initial guess for the variables 
to be solved for; we superpose the respective single-black 
hole Kcrr-Schild quantities, i.e. 



yj = 1, 

n 

m/> = 1+XV^K-l), 



a=l 



a=l 

2/2 
p- r a/ W a ft 1 



(72a) 
(72b) 

(72c) 



where n = 2 and a a and (3 l a are the lapse and shift cor- 
responding to the boosted, spinning Kerr-Schild metrics 
g^j used in the conformal metric <?y . Convergence of the 
elliptic solver and spin are demonstrated in Fig. 10 by 
showing the decreasing constraint violation 8 and differ- 
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FIG. 11: The mass M (upper panel) and dimensionless spin \ 
(lower panel) of one of the holes for Superposed-Kerr-Schild, 
binary-black-hole initial data sets with spins aligned with the 
orbital angular momentum. The mass and spin are plotted 
against Q r [Eq. (40)] for four different choices of the conformal 
spin: S = 0, 0.5, 0.93, and 0.99. Also shown are the data 
sets SKS-0.93-E3— identical to the Q r = 0.28M, S = 0.93M 2 
data set on the solid curve but with lower eccentricity — and 
SKS-Headon; both sets are evolved in Sec. V. The inset in 
the lower panel shows a close-up of the spins as they approach 
unity, with symbols denoting the individual data sets. 



ences in spin with increasing resolution. 

We now turn our attention to the physical properties 
of the SKS initial data sets. Figure 11 shows the horizon 
mass M and the dimensionless spin x of cither black hole 
for the four families of SKS initial data. As expected, we 
find that generally the spin \ increases with increasing 
il r . For each of the SKS-families, we find that the ellip- 
tic solver fails to converge for sufficiently large fl r . We 
suspect that the SKS-families exhibit a turning point, 
similar to the CFMS-single and binary black hole initial 
data shown in Figs. 7 and 9. If this is the case, Fig. 11 
only shows the lower branch of each family, and an ad- 
ditional branch of solutions will be present. Because we 
are satisfied with the spin magnitudes that are possible 
along the lower branch, we do not attempt to find the 
upper branch here. 

In contrast to the CFMS data sets (where the lower 
branch only allowed spins as large as x ^ 0.85), the SKS 
initial data allows spins that are quite close to unity. For 
the different SKS families, we are able to construct initial 
data with spins as large as 

• x ~ 0.95 for SKS-0, 

• x ~ 0.985 for SKS-0.5, 



The constraint violation is y||C||£ 2 + ||C|j^ 2 ||C||^ 2 5ij, where C given by Eq. (73). 



and C l are the residuals of Eqs. (10)-(11) and the L2 norm is 
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FIG. 12: The irreducible mass M\ TT and Euclidean coordinate 
radius r (upper panels) and dimensionless spin \ := S/M 2 and 
spin-extremality parameter £ := S/(2M\ IT ) (lower panels) for 
one of the black holes in the SKS-0.93 (left) and SKS-0.99 
(right) initial-data-set families. These quantities are com- 
puted on two surfaces: i) the apparent horizon (solid lines), 
and ii) the excision boundary of the initial data (dashed lines) . 
Because we enforce that the excision surface is a marginally 
trapped surface, typically the apparent horizon and excision 
boundary coincide. However, if Q r is increased beyond the 
values where \ approaches unity, the apparent horizon lies 
outside of the excision surface. The excision surface can ob- 
tain superextremal spins (C > 1), but only when it is enclosed 
by a subextremal horizon. 

• x ~ 0.998 for SKS-0.93, 

• x ~ 0.9997 for SKS-0.99. 

These spins are far closer to extremal than possible with 
Bowen-York initial data [x < 0.984 (Fig. 2)] or confor- 
mally flat, maximally sliced XCTS initial data [x < 0.85 
or < 0.99 along the lower and upper branch, respectively 
(Fig. 7)]. 

We note that the spins in the SKS binary-black-hole 
initial data families are only weakly dependent on the 
orbital parameters f2o and do- This can be seen from 
the individual data-point labeled SKS-0.93-E3 shown in 
Fig. 11. This data-set uses different values for fio and 
do but is nevertheless close to the family SKS-0.93. The 
initial data sets SKS-0.93-E3 and SKS-HeadOn will be 
discussed in detail in Sec. V. 

The inset of Fig. 11 highlights a remarkable feature of 
the SKS-0.93 and SKS-0.99 families: with increasing fi r , 
the spin initially increases but eventually decreases. Fig- 
ure 12 investigates this behavior in more detail, where 
this effect is more clearly visible in the lower two panels: 
both the spin x and the extremality parameter £ of the 
apparent horizon change direction and begin to decrease. 
For Q r smaller than this critical value, the apparent hori- 
zon finder always converges onto the excision surfaces, 



which by virtue of the boundary condition Eq. (39), are 
guaranteed to be marginally trapped surfaces. As fl r is 
increased through the critical value (at which \ and ( 
change direction), a second marginally trapped surface 
(solid line) splits off from the excision surface (dashed 
line) and moves continuously outward. This can be seen 
in the upper panels of Fig. 12, which plot the minimal 
and maximal coordinate radius and the irreducible mass 
of both the excision surface and the outermost marginally 
trapped surface, which is by definition the apparent hori- 
zon. 

But what about the excision surface? The boundary 
condition Eq. (39) forces the excision surface to be a 
marginally trapped surface, independent of the value of 
il r . For sufficiently large f2 r , however, the excision sur- 
face is surrounded by a larger marginally trapped surface 
and thus is not the apparent horizon. The dashed lines 
in Fig. 12 present data for the excision surface. These 
lines continue smoothly across the point where the sec- 
ond marginally trapped surface forms. The extremal- 
ity parameter £ for the excision surface continues to in- 
crease and eventually becomes larger than unity; the ex- 
cision surface can then be thought of as having a superex- 
tremal spin. However, for the outer marginally trapped 
surface — the true apparent horizon — the extremality pa- 
rameter always satisfies £ < 1 . The irreducible mass Afj rr 
of this surface increases faster than the spin, and there- 
fore C = S/(2M? T ) decreases with increasing Q r . 

One might interpret these results as support of the 
cosmic censorship conjecture. The XCTS boundary con- 
ditions (39) and (40) control the location and the spin of 
the excision surface. By appropriate choices for the shift 
boundary condition (40), we can force the excision sur- 
face to become superextremal. However, before this can 
happen, a new horizon appears, surrounding the excision 
surface and hiding it from "our" asymptotically flat end 
of the spacetime. The newly formed outer horizon always 
remains subextremal. 



C. Suitability for evolutions 

In the previous sections, we have constructed a wide 
variety of binary-black- hole initial data sets. To get some 
indication about how suitable these are for evolutions, 
we consider the initial time-derivatives of these data sets, 
dt9ij and dtKij. Recall that solutions of the XCTS equa- 
tions give a preferred initial lapse and shift for the evolu- 
tion of the initial data; hence, the time derivatives dtgij 
and dtKij can be computed by simply substituting the 
initial data into the ADM evolution equations. We ex- 
pect initial data with smaller time-derivatives to be closer 
to quasi-equilibrium and to have less initial spurious ra- 
diation. 

Figure 13 presents the L2 norms of the time deriva- 
tives, ||<9t<?ij llj^ and ||dtiQj|| i2 where the L2 norm of a 
tensor Tijk...(x) evaluated at N gridpoints xi is defined 



17 



— CFMS 
*-* SKS-0.0 
x-x SKS-0.5 
->-+ SKS-0.93 
SKS-0.99 
SKS-0.93-E3 
SKS-Headon 





FIG. 13: The time derivatives of the metric (left panel) and 
extrinsic curvature (right panel). In the superposed-Kerr- 
Schild (SKS) data sets, ||3t.KVj || i2 has minima near values of 
Q r for which the dimensionless spin \ i s approximately equal 
to the spin S of the conformal metric (cf. Fig. 11). On the up- 
per branch of the conformally-flat, maximally-sliced (CFMS) 
excision data, where the spin is \ > 0.83 (Fig. 9), the time 
derivatives become much larger than the SKS time deriva- 
tives. The data sets SKS-0.93-E3 (with X ~ S = 0.93) and 
SKS-Headon (with x ~ S = 0.97) are evolved in Sees. VC- 
VD; the time derivatives are significantly lower for the set 
SKS-Headon because of the larger coordinate separation of 
the holes (d = 100 vs. d = 32). 



as 



\T, 



i=0 



where 



T 



-Mjfe-" -me' ' j'k' 



_S iif 5 JJ, S kk ' , 



(73) 



(74) 



Figure 13 shows that generally dtKij is larger than 
dtgtj ■ This has also been found in previous work, e.g. [65] , 
and is not surprising, because the XCTS formalism al- 
lows some control over the time derivative of the metric 



through the free data Ui 



whereas there is less 



control of dtKij. We note that for CFMS data, the time 
derivatives are larger and grow more rapidly with x than 
for SKS data; in particular, the time derivatives on the 
upper branch arc ~ 10 times larger than for SKS-initial 
data, suggesting that these data are much farther from 
equilibrium. 

In the SKS case, the time derivatives of Kij have lo- 
cal minima at particular values of Q r ; comparison with 
Fig. 11 gives spins \ at these minima of ||c?iif< 
follows: 

• SKS-0.5: Q r « 0.1, X ~ 0.45, 

• SKS-0.93: fl r w 0.28, X ~ 0.93, 

• SKS-0.99: fl r w 0.34, X « 0.98. 



»J \\L2 



a.s 



Note that these minima occur at values of fl r such that 
X ~ S/M 2 ; that is, transients in the initial data and 
presumably the spurious radiation are minimized when 
the conformal spin and AKV spin are consistent. For 
this reason, we conclude that SKS initial data with x ~ 
S/M 2 is preferable; this is the type of initial data we will 
evolve in the next section. 

Also note that minimizing the spurious radiation has 
purely numerical advantages: the spurious radiation typ- 
ically has finer structure (and thus requires higher reso- 
lution) that the physical radiation. If such radiation is 
minimized, the numerical evolutions may require less res- 
olution and will be more efficient. Conformally-curved 
initial data has been found to reduce the amount of spu- 
rious radiation in Refs. [39, 66]. 



V. EXPLORATORY EVOLUTIONS OF 
SUPERPOSED KERR-SCHILD (SKS) 
INITIAL DATA 



So far, we have confined our discussion to black hole 
spins in the initial data. In this section, we compare the 
initial spin to the value to which the spin relaxes after the 
initial burst of spurious radiation, when the holes have 
settled down. Recall, for instance, that for Bowen-York 
puncture initial data with spins close to the maximal 
possible value [x(t = 0) « 0.98], the spins quickly re- 
lax by about A% w 0.05 to a maximal possible relaxed 
value of xOroiax) ~ 0.93 (cf. [25]). While the SKS data 
presented in Sec. IV B can achieve larger initial spins 
[x(t — 0) = 0.9997] than conformally-flat puncture data, 
only evolutions can determine Ax and x(^reiax)- 

Therefore, in this section we perform brief, exploratory 
evolutions of some SKS initial data sets to determine Ax 
for those data sets. 9 Besides determination of x(^reiax), 
these evolutions will also allow us to demonstrate that the 
technique of eccentricity reduction developed in Rcf. [49] 
is applicable to SKS initial data as well as to compare 
the spin measures defined in Appendices A and B. The 
focus here lies on initial data, and we evolve only long 
enough for our purposes. Longer simulations that con- 
tinue through merger and ringdown are the subject of 
ongoing research. 

This section is organized as follows. In Sec. VA, 
we summarize the evolution code that we will use. In 
Sec. VB, we perform eccentricity reduction on one of 
the data sets in the SKS-0.93 family, which corresponds 
to an orbiting binary black hole with equal masses and 
equal spins (of magnitude x ~ 0.93) aligned with the or- 
bital angular momentum. Then, in Sec. VC, we evolve 
the resulting low-eccentricity data set (labeled SKS-0.93- 



9 Note that there is no universal value of Ax — it will differ for 
different initial data sets, even within the same family of initial 
data. 
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E3). Finally, In Sec. VD, we evolve a head-on plunge of 
SKS initial data (labeled SKS-Headon) representing two 
widely-separated black holes with initial spins of magni- 
tude x = 0.970 and direction normal to the equatorial 
plane. 



A. Description of evolution code 



The initial data are evolved using the Caltech-Cornell 
pscudospectral evolution code SpEC [50]. The details of 
the evolution methods, equations, and boundary condi- 
tions that we use are the same as those described in 
Ref. [67]. The singularities are excised, with the exci- 
sion surfaces chosen to lie slightly inside the black hole 
horizons. Note that whereas Ref. [67] excises coordinate 
spheres inside the black holes' apparent horizons, here 
we use Lorentz-contracted ellipsoidal excision boundaries 
which are adapted to the shape of the initial apparent 
horizons. 

The highest-resolution initial data set (with N as 85 3 
gridpoints) is interpolated onto evolution grids labelled 
Nl, N2, and A3 with approximately 61 3 , 67 3 , and 74 3 
gridpoints, respectively. The outer boundary is at a co- 
ordinate radius of r = 32d for the orbiting simulation 
discussed in Sees. VB and VC and at r = 14d for the 
head-on simulations discussed in Sec. VD. This trans- 
lates to about r = 450i?ADM and r = 620-Eadm for the 
orbiting and head-on simulations, respectively. As in ear- 
lier simulations [49, 50, 67], a small region of the evolu- 
tion grid lies inside the horizon and is not covered by the 
initial data grid; we extrapolate ip, atp, and (¥" into this 
region and then compute gij and Kij . 



B. Eccentricity removal for orbiting SKS-binaries 

We obtain initial data with small orbital eccentric- 
ity using the iterative method of Ref. [49] , as refined in 
Ref. [67], applied here for the first time to binary-black- 
hole data with rapid spin. In this method, the choice 
of flo and do for the next iteration are made so that if 
the orbit were Newtonian, the eccentricity would vanish. 
For the non-Newtonian orbit here, successive iterations 
succeed in reducing the orbital eccentricity. 

This procedure is based on the proper separation s 
between the apparent horizons, measured along a coordi- 
nate line connecting the geometric centers of the apparent 
horizons. The time derivative ds/dt is fitted to a five- 
parameter curve that, together with the initial proper 
separation s(t = 0) is used to define the eccentricity e 
and to define improved values for SIq and do- Specifi- 
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FIG. 14: Color online. Eccentricity reduction for evolutions 
of superposed-Kerr-Schild binary-black-hole initial data. The 
proper separation s (upper panel) and its time derivative 
ds/dt (lower panel) are plotted for initial data sets SKS-0.93- 
E0, -El, -E2, and -E3, which have successively smaller ec- 
centricities e. All evolutions are performed at resolution Nl. 
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Heuristically, the eccentricity is embodied by the oscil- 
lating part of ds/dt. 

Figure 14 illustrates the eccentricity reduction for one 
of the data sets in family SKS-0.93. Plotted are the 
proper separation s and its derivative ds/dt for evolu- 
tions of several initial data sets (summarized in Table I): 

• set SKS-0.93-E0, which is identical to the set in 
family SKS-0.93 with fl r = 0.28 (Fig. 11); 

• set SKS-0.93-E1, which is the same as SKS-0.93- 
E0 except that the orbital frequency f2o is manually 
adjusted to lower the orbital eccentricity somewhat; 
and 

• sets SKS-0.93-E2 and SKS-0.93-E3, which arc suc- 
cessive iterations (starting from set SKS-0.93-E1) 
of the eccentricity-reduction scheme Eqs. (75). 

The ad hoc adjustment of f^o was somewhat effective, re- 
ducing e by about 50%. The subsequent iterations using 
Eqs. (75) reduced e by factors of about 5 and 8, respec- 
tively. Surprisingly, the lowest eccentricity, correspond- 
ing to a smooth inspiral trajectory is obtained with a 
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FIG. 15: Convergence test of the evolution of the initial data 
set SKS-0.93-E3. Shown are evolutions on three different res- 
olutions, Nl, N2, and 7V3, with 7V3 being the highest reso- 
lution. The top panel shows the approximate-Killing-vector 
(AKV) spin of one of the holes as a function of time, with 
the top inset showing the spin's initial relaxation; the bottom 
panel shows the constraint violation as a function of time. 
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FIG. 16: A comparison of different definitions of the spin. 
The top panel shows the spin as a function of time for sev- 
eral different measures of the spin; the bottom panel shows 
the fractional difference between xakv and alternative spin 
definitions. Note that for t < 30-Eadm, the time-axis has a 
different scaling to make the initial transients visible. 



positive do = 3.332 x 10~ 4 . This is not due to insuffi- 
cient resolution; for SKS-0.93-E3, we have verified that 
we obtain the same eccentricity e ~ 0.001 for all three 
numerical resolutions Nl, N2, N3. 

Note that we choose to stop the evolutions at about t = 
670i?ADM, which corresponds to about 1.9 orbits; this is 
sufficient for reducing the eccentricity and for measuring 
Ax- In the next subsection, we discuss the evolution of 
the low-eccentricity set SKS-0.93-E3 in detail, focusing 
on the relaxation of the spin 



C. Low-eccentricity inspiral with \ ~ 0.93 

We evolved the data-set SKS-0.93-E3 at three different 
numerical resolutions for a duration of about 670-Eadm, 
corresponding to about 1.9 orbits. From post-Newtonian 
theory [68], we estimate that this simulation would pro- 
ceed through about 20 orbits to merger. 

Figure 15 presents a convergence test for this run. The 
lower panel of Fig. 15 shows the normalized constraint 
violation (see Eq. (71) of Ref. [69] for the precise defini- 
tion.) While the constraints are small, the convergence 
seems poor until t sw 500i?ADM- For this time-period the 
constraint violations at high resolution N3 are dominated 
by the outgoing pulse of spurious radiation — i.e. far away 
from the black holes — which we have not attempted to 
adequately resolve. At t ~ 500-EadM; the pulse of spuri- 
ous radiation leaves the computational domain through 
the outer boundary; afterwards, the constraints decrease 
exponentially with increasing resolution, as expected. 

The upper panel of Figure 15 shows the AKV spin 
Xakv = S/M 2 for the three runs with different reso- 
lutions Nl, N2, and N3. Based on the difference be- 



tween N2 and N3, the spin of the evolution N3 should 
be accurate to a few parts in 10 . For the time-interval 
5 < tj Eaum < 670, the measured spin on resolution iV3 
is consistent with begin constant within its estimated ac- 
curacy. Very early in the simulation, t < 5-Eadm, the 
spin x changes convergently resolved from its initial value 
X(t = 0) = 0.927 48 to a relaxed value xOreiax) = 0.927 14 
(see inset of Fig. 15). Therefore, for SKS-0.93-E3, we find 
Ax = 0.000 34. 

Contrast this result with the evolution of a binary 
black hole puncture initial data set with large spins, 
which is reported in Ref. [25]: for that particular evo- 
lution, x (t = 0) = 0.967, x(irdax) = 0.924, i.e. A x = 
0.043, more than a factor 100 larger than for the evolu- 
tion of SKS-0.93-E3 reported here. This comparison is 
somewhat biased against the puncture evolution in [25], 
which starts at a smaller separation possibly resulting in 
larger initial transients. However, even in the limit that 
the black holes are infinitely separated (i.e., in the singlc- 
black-hole limit), the spins in Bowen-York puncture data 
relax to values near £j = Jadm / E\ DM \ to achieve a fi- 
nal spin of x(^reiax) ~ 0.93, the initial spin of Bowen-York 
data must be x(t = 0) w 0.98 (cf. Fig. 2 of Ref. [25]). We 
conclude that the spin relaxes by a much smaller amount 
in the SKS case than in Bowen-York puncture or inver- 
sion symmetric data. 

Figure 15 and the discussion in the previous paragraph 
only addresses the behavior of the AKV spin, where the 
approximate Killing vectors are computed from the min- 
imization problem [cf. Eq. (A10)]. We now compare 
the different spin-definitions we present in Appendices A 
and B. Figure 16 compares these different definitions of 
the black hole spin for the N3 evolution of initial data set 
SKS-0.93-E3. Shown are the AKV spin of one hole in the 



20 



binary, the scalar curvature (SC) spins XscT an d Xsc* °f 
Appendix B [Eqs. (B2a) and (B2b)], and also the spin ob- 
tained by using Eq. Al with a coordinate rotation vector 
instead of an approximate Killing vector (which we call 
the "coordinate spin" here). After the holes have relaxed, 
the SC spins track the AKV spin more closely than does 
the coordinate spin. However, during very early times, as 
the holes are relaxing and the horizon shape is very dis- 
torted, the SC spins show much larger variations. Con- 
sequently, the SC spin is a poorer measure of the spin at 
early times than even the coordinate spin. 



D. Head-on plunge with % ~ 0.97 

In the previous subsection, we have seen that for SKS 
binary-black hole-initial data with x = 0.93, the initial 
spins change by only a few parts in 10 4 . A spin % « 0.93 
is roughly the largest possible equilibrium spin that is 
obtainable using standard conformally-flat, Bowcn-York 
puncture data (cf. the discussion at the beginning of 
Sec. V). We now begin to explore binary-black- hole sim- 
ulations with spin-magnitudes that are not obtainable 
with Bowen-York initial data methods. 

We construct and evolve SKS binary-black-hole data 
for a head-on plunge of two equal mass black holes with 
spins of equal magnitude x — an d with the spins 
orthogonal to the line connecting the black holes. This 
data set, labelled SKS-Headon, is summarized in Table I 
and was briefly discussed in Sec. IV B, cf. Figs. 10, 11 
and 13. As for the orbiting evolution SKS-0.93-E3, we 
adjust the rotation parameter Q r so that conformal spin 
S/M 2 and AKV spin x are approximately equal. Start- 
ing such a simulation at close separation results in rapid 
coordinate motion of the apparent horizons during the 
first few £adm of the evolution. These motions are cur- 
rently difficult to track with our excision code; therefore, 
we begin at a larger separation d than we used in the 
nearly-circular data sets described previously. 

Figure 17 presents a convergence test of the constraints 
(lower panel) and the AKV spin xakv (upper panel) 
during the subsequent evolution. Again, we are inter- 
ested in the initial relaxation of the spins; therefore, we 
choose to stop evolution at t rts 120-Eadm- During this 
time, the black hole proper separation decreased from 
s(t = 0) = 47.6£ A dm to s(t = 120) = 44.1£ A dm- 

During the first ~ 10-Eadm, Xakv shows (a numeri- 
cally resolved) decrease of about 3 x 10 -5 ; this change 
arises due to initial transients as the black holes and the 
full geometry of the spacetime relax into an equilibrium 
configuration. Subsequently, the spin remains constant 
to within about 10 -4 , where these variations are domi- 
nated by numerical truncation error. 

Figure 18 compares our various spin-measures for the 
head-on simulation. Interestingly, the spin Xcoord com- 
puted from coordinate rotation vectors agrees much bet- 
ter with xakv than for the SKS-0.93-E3 evolution, per- 
haps because the black holes here are initially at rest. 
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FIG. 17: Convergence test of the head-on evolution SKS- 
HeadOn. Shown are evolutions at three different resolutions, 
Nl, N2, and N3, with iV3 being the highest-resolution. The 
top panel shows the approximate-Killing-vector (AKV) spin 
of one of the holes as a function of time; the bottom panel 
shows the constraint violations as a function of time. 
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FIG. 18: A comparison of various measures of the spin for 
the head-on evolution of data set SKS-Headon, which is a 
plunge of two equal-mass black holes with with parallel spins 
of magnitude xakv = S/M 2 = 0.970 pointed normal to the 
equatorial plane. The top panel shows various measures of the 
spin as a function of time, and the bottom panel shows the 
fractional difference between the approximate-Killing-vector 
(AKV) spin xakv and alternative spin definitions. 



The scalar-curvature (SC) spins Xsc" an d XscfS de- 
rived from the scalar curvature of the apparent horizon 
[Eqs. (B2a) and (B2b)], show some oscillations at early 
times; after the initial relaxation, the SC spin agrees with 
the AKV spin to about 1 part in 10 4 . 
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VI. DISCUSSION 
A. Maximal possible spin 

In this paper, we have examined a variety of meth- 
ods for constructing black hole initial data with a par- 
ticular emphasis on the ability to construct black holes 
with nearly-extremal spins. These arc spins for which 
the dimensionless spin \ = S/M 2 and spin-extremality 
parameter ( = S/ (2M? r ) are close to unity. 

When discussing black hole spin, one needs to distin- 
guish between the initial black hole spin and the relaxed 
spin of the holes after they have settled down. Using 
conformally-fiat Bowen-York (BY) data (both puncture 
data or inversion symmetric data) for single black holes, 
the largest obtainable spins are \ sa 0.984, £ sa 0.833 (cf. 
Ref. [63] and Fig. 2). With conformally-flat, maximally- 
sliced (CFMS), quasi-equilibrium cxtendcd-conformal- 
thin-sandwich (QE-XCTS) data, we are able to obtain 
initial spins as large as x ~ 0.99, C ~ 0.87 for single 
black holes (Fig. 7). The limitations of BY puncture data 
and CFMS QE-XCTS data are already present when con- 
structing highly spinning single black holes; therefore, we 
expect the methods to be able to construct binary-black- 
hole data with similar spins as for single holes — i.e., up 
to about 0.98. Construction of CFMS QE-XCTS binary- 
black-hole initial data confirms this conjecture (compare 
Fig. 9 with Fig. 7). 

For superposed-Kerr-Schild (SKS) initial data, the sit- 
uation is different. For single black holes, SKS data re- 
duce to the analytical Kerr solution, without any limi- 
tations on the spin magnitude. Thus limitations of SKS 
data will only be visible for binary-black hole configura- 
tions. As Sections IV and V show, however, those limita- 
tions are quite minor. SKS data can indeed achieve ini- 
tial spins that are much closer to extremality than what 
is possible with BY data or CFMS QE-XCTS data; we 
have explicitly demonstrated this by constructing SKS 
data for binary black holes with \ ~ 0.9997, ( ~ 0.98, as 
can be seen from Figs. 11 and 12. 

As the black hole spacetimes settle into equilibrium 
and emit spurious gravitational radiation, the initial spin 
X decreases to a smaller relaxed spin x(^rciax)- Thus 
an interesting quality factor for high-spin black hole ini- 
tial data is Ax = x(t = 0) — x(^reiax) [Eq. (3)] consid- 
ered as a function of the relaxed spin. The magnitude 
of Ax is indicative of the amplitude of any initial tran- 
sients, whereas the maximally achievable x(^-eiax) gives 
the largest possible spin which can be evolved with such 
initial data. Figure 19 presents this plot, with the circle 
and cross representing the two evolutions of SKS data 
which were described in Sec. V. 

We have not evolved high-spin puncture data, nor 
high-spin CFMS-XCTS data; therefore, we do not know 
precisely Ax for these initial data. We estimate Ax for 
puncture data by noting that evolutions of single-hole, 
BY puncture data with large spins show [25] that the 
black hole spin x '■= S/M 2 relaxes approximately to the 
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FIG. 19: The change in black hole spin \ during the initial 
relaxation of black hole initial data plotted as a function of 
the black-hole spin after relaxation. The SKS initial data 
constructed in this paper have smaller transients and allow 
for larger relaxed spins. 

initial value of Ej := Jadm/^adm- Therefore, for BY 
puncture data, we approximate 

A x ^ej- X (t = 0), (76) 

X(*relax) « Bj. (77) 

This curve is plotted in Fig. 19. Because high-spin singlc- 
black-hole, CFMS QE-XCTS initial data and BY punc- 
ture data have quite similar values of x{t = 0) and ej, as 
well as similar embedding diagrams (cf. Fig. 8), we con- 
jecture that Eqs. (76)-(77) are also applicable to CFMS 
QE-XCTS data. This estimate is also included in Fig. 19. 
We see that both types of initial data result in a Ax of 
similar magnitude which grows rapidly with Xrciaxed- 

Perhaps the most remarkable result of Fig. 19 is the 
extremely small change in black hole spin during the re- 
laxation of SKS initial data, even at spins as large as 
X = 0.97. The small values of Ax combined with the 
ability to construct initial data with initial spins x(t = 0) 
as large as 0.9997 (cf. Fig. 11) makes it highly likely that 
SKS initial data are capable of constructing binary black 
holes with relaxed spins significantly closer to unity than 
0.97. Evolutions of initial data with spins x much closer 
to unity, i.e., farther into the regime that is inaccessi- 
ble to conformally-flat data, are a subject of our ongoing 
research. 

In summary, the two main results of this paper are as 
follows: 

• SKS initial data can make binary black holes that 
initially have nearly-extremal spins, and 

• for SKS initial data, the relaxed spin is quite close 
to the initial spin, even when the spin is large. 
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B. Additional results 

While working toward the main results discussed in 
the previous subsection, we have also established several 
additional interesting results. We have considered spin- 
ning, single-black-hole, puncture data which is identical 
to single-black-hole, spinning, inversion-symmetric data. 
Using this correspondence and our accurate spectral el- 
liptic solver, we revisited the relation between black-hole 
spin x, specific total angular momentum of the space- 
time Sj, and the spin-parameter S for BY puncture data, 
and established in Fig. 2 that both \ an d £j approach 
their limits for S — > oo as power-laws, cf. Eqs. (49) 
and (50). We have also extended the analytical analysis 
of Dain, Lousto, and Zlochowcr [25] of the throat region 
of high-spin puncture data toward more quantitative re- 
sults, including the precise amplitudes of the conformal 
factor, throat circumference and throat length, as well as 
their scaling with spin-parameter S and puncture mass 
m p [Eqs. (54)-(58)]. Furthermore, Ref. [25] implicitly 
assumed that the throat-region is approximately spher- 
ically symmetric; our Fig. 5 presents explicit evidence 
in support of this assumption, but also shows that the 
throat is not precisely spherically symmetric. 

We have also examined high-spin QE-XCTS initial 
data employing the common approximations of confor- 
mal flatness and maximal slicing (CFMS). With increas- 
ing angular frequency Q r of the horizon, we discover non- 
unique solutions. Thus, the non-uniqueness of the XCTS 
equations can not only be triggered by volume terms (as 
in [34]) but also through boundary conditions [in this 
case, by Eq. (40)]. Interestingly, CFMS QE-XCTS data 
appears to be very similar to BY puncture data, in regard 
to nearly-extremal spins. Both data formalisms result in 
similar maximal values of x(t = 0) and sj (Figs. 2, 7 
and 19) and have embedding diagrams which develop a 
lengthening throat as the spin is increased (Fig. 8). 

We also have found an interesting property of the hori- 
zon geometries for SKS data, which one might interpret 
as support of the cosmic censorship conjecture. Specif- 
ically, we find that by increasing Cl r sufficiently, we can 
in fact force the excision boundaries of the initial data 
to be "horizons" (i.e. marginally trapped surfaces) with 
superextremal spin (£ > 1). However, these superex- 
tremal surfaces are always enclosed by a larger, subex- 
tremal (£ < 1) apparent horizon. 

To measure black hole spins, we have employed and 
compared several different techniques to measure black 
hole spin. Primarily, we use a quasilocal spin definition 
based on (approximate) Killing vectors [Eq. (Al)]. This 
formula requires the choice of an "approximate" Killing 
vector, and we have used both straightforward coordinate 
rotations to obtain Xcoord and solved Killing's equation 
in a least-squares sense to obtain xakv (see Appendix A 
or details). Furthermore, we introduced a new technique 
to define black-hole spin which docs not require choice 
of an approximate Killing vector and is invariant under 
spatial coordinate transformations and transformations 



associated with the boost gauge ambiguity of the dynam- 
ical horizon formalism. This new technique is based on 
the extrema of the scalar curvature of the apparent hori- 
zon. Figures 16 and 18 show that all four spin measures 
agree to good precision, but differences are noticeable. 
The spin-measures based on the horizon curvature ex- 
hibit more pronounced variations during the initial tran- 
sients, and the quasilocal spin based on coordinate rota- 
tions is off by several tenths of a percent. The quasilocal 
spin based on approximate Killing vectors xakv has the 
smallest initial variations. 

Finally, we would like to point out that a modified 
version of the SKS-initial data has been very successfully 
used to construct black hole- neutron star initial data [70] . 
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APPENDIX A: QUASILOCAL SPIN USING 
APPROXIMATE KILLING VECTORS (AKV 
SPIN) 

In this appendix and the one that follows, we address 
the task of defining the spin of a dynamical black hole, 
given gij and Kij. We use two different measures. The 
first, defined here, is a standard quasilocal angular mo- 
mentum defined with approximate Killing vectors which 
correspond to approximate symmetries of a black hole's 
horizon. The second measure, defined in Appendix B, 
infers the spin from geometrical properties (specifically, 
from the intrinsic scalar curvature R) of the apparent 
horizon, assuming that the horizon is that of a single 
black hole in equilibrium, (i.e., that the horizon is that 
of a Kerr black hole). Note that quantities relating to 
the geometry of the two-dimensional apparent horizon 
surface H. are denoted with a ring above them, to avoid 
confusion with the analogous quantities on the spatial 
slice, S. 
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It has become standard in the numerical relativity 
community to compute the spin angular momentum of 
a black hole with the formula [71-73] 

S= ~kf ^ sJKlJ dA ' (A1) 

where s l is the outgoing normal of TL embedded in £ 
and </> is an "azimuthal'' vector field, tangent to Ti.. The 
azimuthal vector field </> carries information about the 
"axis" about which the spin is being computed. There 
are, however, far more vector fields on a two sphere than 
there are axes in conventional Euclidean space. We must 
find suitable criteria for fixing these azimuthal vector 
fields in numerical simulations, so that they reduce to 
the standard rotation generators when considered on a 
metric sphere. 

Because angular momentum is generally thought of as 
a conserved charge associated with rotation symmetry — 
and indeed the quantity given in (Al) can be shown to 
be conserved under time evolution [71, 73] when is 
a Killing vector of the dynamical horizon worldtube — it 
makes sense to consider Killing's equation to be the es- 
sential feature of the azimuthal vector field. If a Killing 
vector on a dynamical horizon is tangent to each (two- 
dimensional) apparent horizon, then the vector field must 
be a Killing vector of each apparent horizon. However 
in a general spacetime, on an arbitrary apparent hori- 
zon, there is no reason to expect any Killing vectors to 
exist. So in the cases of most interest to numerical rel- 
ativity, when there are no true rotation symmetries, we 
must relax the symmetry condition and find those vector 
fields that come "closest" to generating a symmetry of 
the apparent horizon. In other words, we seek optimal 
"approximate Killing vectors" of the apparent horizon. 

In [74], a practical method for computing approxi- 
mate Killing vectors was introduced, which has since 
been applied on numerous occasions, e.g. [18, 19, 29, 75] 
This method involves integrating the Killing transport 
equations along a predetermined network of coordinate 
paths. The resulting vector field is guaranteed to be a 
Killing vector field if such a field exists and coincides 
with the computed field at any point on the network. 
However if no true Killing field exists, the integral of 
the Killing transport equations becomes path dependent. 
This means that the computed vector field will depend 
in an essential way on the network of paths chosen for 
the integral. Perhaps even more serious, if there is no 
true Killing field, then the transport of a vector around 
a closed path will not necessarily be an identity map. As 
a result, the computed vector field cannot be expected 
to reduce to any smooth vector field in the limit that 
the network becomes more refined. This kind of approx- 
imate Killing vector field is simply not mathematically 
well-defined in the continuum limit. 

Here we will describe a kind of approximate Killing 
vector field that, as well as having a well-defined contin- 
uum limit, is actually easier to construct than those of 



the Killing transport method, at least in our particular 
code. Our method is extremely similar to that described 
by Cook and Whiting [41], but was actually developed 
independently by one of the current authors [76] . 

1. Zero expansion, minimal shear 

Killing's equation, 

D {A cj) B) = 0, (A2) 

has two independent parts: the condition that 4> be 
expansion- free, 

6 := g AB D A <t> B = 0, (A3) 
and the condition that it be shear-free, 

°ab ■= D {A (f> B) - -g AB @ = 0, (A4) 

where uppercase latin letters index the tangent bundle to 
the two-dimensional surface, g AB is the metric on that 
surface, and D A is the torsion-free covariant derivative 
compatible with that metric. 

When constructing approximate Killing vectors, a 
question arises: which condition is more important, zero 
expansion or zero shear? Shear-free vector fields (confor- 
mal Killing vectors) arc simply coordinate rotation gener- 
ators in the common case of coordinate spheres in a con- 
formally flat space. They are therefore readily available 
in that context. A very interesting and systematic ap- 
proach to their use has been given by Korzynski [77] , and 
they have been used in the construction of conformally 
flat binary black hole initial data sets [28, 29]. However, 
in the case of a general surface in a general spatial slice, 
the conformal Killing vectors are not known a priori, and 
they are more difficult to construct than expansion-free 
vector fields. Expansion-free vector fields have the addi- 
tional benefit of providing a gauge-invariant spin measure 
on a dynamical horizon [73] 10 , so we restrict attention to 
the expansion-free case. 

Any smooth, expansion-free vector field tangent to a 
topological two-sphere can be written as 

4> A = e AB D B z, (A5) 

where e AB is the Levi-Civita tensor and z is some smooth 
potential function. 



The dynamical horizon is essentially the world tube foliated by 
the apparent horizons. The gauge freedom is that of extending 
the foliation off of this world tube. The fact that this gauge 
invariance occurs when (f> is expansion-free can most easily be 
shown by expressing the factor s J Kij in Eq. (Al) in terms of the 
ingoing and outgoing null normals to the two-surface. The boost 
freedom in these null normals has no effect on the spin when <j> 
is expansion-free. 
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We assume that the function z has one local maximum, 
one local minimum, and no other critical points. This 
is equivalent to the assumption that the orbits of 4> are 
simple closed loops. In order for <p A (j)A to have the proper 
dimensions, z must have dimensions of area. For the case 
of the standard rotation generators of the metric two- 
sphere, the three z functions are the three £ = 1 spherical 
harmonics, multiplied by the square of the areal radius 
of the sphere. 

Within this space of expansion-free vector fields, we 
would now like to minimize the following positive-definite 
norm of the shear: 

lk|| 2 := I a BC a BC dA. (A6) 

Substituting Eq. (A5) for (j> in this expression and inte- 
grating twice by parts, ||cx|| 2 takes the form of an expec- 
tation value: 

||cr|| 2 = I zHzdA, (A7) 

where H is the self-adjoint fourth-order differential oper- 
ator defined by 

Hz = D 4 z + RD 2 z + D A R D A z, (A8) 

and D 2 is the Laplacian on the (not necessarily round) 
sphere, is its square, and R is the Ricci scalar cur- 
vature of the sphere. In our sign convention, R = 2 on 
the unit sphere, so we can immediately see that Hz = 
when z is an £ — 1 spherical harmonic, and therefore that 
their associated vector fields are shear-free. 

It is now tempting to minimize the functional ||er|| 2 in 
(A7) with respect to z. However, doing so will simply 
return the condition that z lie in the kernel of H. If 
there are no true Killing vectors, this will mean that z 
is a constant, and therefore that <f> vanishes. We need to 
restrict the minimization procedure to cases that satisfy 
some normalization condition. In this case, we require 
that the norm of the vector field, 

I cf) A <f) A dA, (A9) 

take some given positive value. This restriction can be 
made with the use of a Lagrange multiplier. Specifically, 
the functional we wish to minimize is 

I[z] := j> zHz dA + \(^j> D A z D A z dA - Nj (A10) 

for some yet undetermined positive parameter N . Note 
that A is the Lagrange multiplier and we have made use 
of the fact that Eq. (A5) implies that (j> ■ <j> — Dz ■ Dz. 
Minimizing the functional / with respect to z returns a 
generalized eigenvalue problem: 

Hz = \D 2 z. (All) 



It is at this point that we can most easily clarify the dif- 
ference between our construction of approximate Killing 
vectors and that of Cook and Whiting in [41]. The differ- 
ence lies in the choice of norm in which the minimization 
problem is restricted. Rather than fixing the norm (A9) 
to take some fixed value in the minimization, Cook and 
Whiting instead fix the dimcnsionlcss norm: 

(f Rcj) A <f)A dA. (A12) 

In general, we see no particular reason to prefer either 
norm over the other, but for the current purposes we 
have at least an aesthetic preference for (A9), which is 
positive-definite even at high spin, whereas (A12) is not, 
because the scalar curvature R of the horizon becomes 
negative near the poles at high spin. If the norm (A9) in 
Eq. (A10) is replaced by (A12), the result is the problem 
described in [41]: 

Hz = X(RD 2 z + D A R D A z). (A13) 

In our numerical code, we discretize (All) (or, option- 
ally, (A13), but not for any results published here) and 
solve the resulting linear algebra problem with a LA- 
PACK routine [78]. Note, however, one technical pecu- 
liarity: the operators H and D 2 in (All) share a kernel, 
the space of constant functions. This means that this 
generalized eigenvalue problem is singular, a fact that 
can cause considerable difficulties for the numerical so- 
lution [79]. The same can be said of (A13). For our 
purposes, this complication is easily evaded. Since we are 
working with a spectral code, it is easiest to discretize the 
problem using expansion into the spectral basis functions 
(coordinate spherical harmonics). When this is done, the 
space of constant functions — the shared kernel of the two 
operators — is simply the span of a single basis function: 
the constant, Yqq. This basis function can easily be left 
out of the spectral expansion, and thereby removed from 
the numerical problem. 

Expansion into coordinate spherical harmonics has an- 
other practical advantage. As noted earlier, for metric 
spheres in standard coordinates, the Killing vectors arise 
when z is given by an £ = 1 spherical harmonic. Thus, 
assuming our horizon is nearly round, and noticeably so 
in the given coordinates, the lowest basis functions (the 
( = 1 spherical harmonics) should nicely approximate 
the intended eigenfunctions. The higher basis functions 
should simply provide small corrections. 

In summary, the approach that we take to finding ap- 
proximate Killing vectors begins with a spectral decom- 
position of Eq. (All). This problem, of course, provides 
as many eigenvectors as there are elements of the spec- 
tral decomposition. We restrict attention to the three 
eigenvectors with smallest eigenvalues (ignoring the vec- 
tor corresponding to the constant eigenfunction, which 
is physically irrelevant and removed from discretization) , 
as these are the ones corresponding to vector fields with 
the smallest shear, and at least for spheres that are only 
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slightly deformed, the orbits of these vector fields are 
smooth closed loops. 

It must be noted that only the eigenvector with the 
smallest eigenvalue corresponds to a vector field with 
strictly minimum shear: even locally, all other eigen- 
vectors are saddle points of the minimization problem. 
The three of them taken together, however, provide a 
geometrically-defined subspace of the vector space of 
expansion-free vector fields, a natural generalization of 
the rotation generators on metric spheres. Using these 
three vector fields (normalized as described in the next 
subsection), one can define "components" of the spin 
angular momentum of a black hole 11 , and from these 
components infer the spin around an arbitrary axis or 
even a spin "magnitude" using a metric on this three- 
dimensional space of generalized rotation generators. In 
practice, we have found no need to go quite so far. As 
mentioned in [41], the approximate Killing vectors gen- 
erally adapt themselves so well to the horizon that one of 
the components is much larger than the other two, so this 
is considered the spin magnitude, and the associated ap- 
proximate Killing vector is considered to define the spin 



2. Normalization 

Solutions to the cigcnproblem (All) can only deter- 
mine the approximate Killing vectors up to a constant 
scaling. Fixing this scaling is equivalent to fixing the 
value of N in (A10). The standard rotation generators 
of metric spheres are normalized such that, when con- 
sidered as differential operators along their various or- 
bits, they differentiate with respect to a parameter that 
changes by a value of 2ir around each orbit. Naively 
one would like to fix the normalization of approximate 
Killing vectors in the same way, but a subtlety arises: 
we can only rescale the vector field by a fixed, constant 
value. Rescaling differently along different orbits would 
introduce extraneous shear and would remove the vec- 
tor field from the pure eigenspace of (All) in which it 
initially resided. If an approximate Killing vector field 
has different parameter circumferences around different 
orbits, then it is impossible to rescale it such that the pa- 
rameter distance is 2ir around every orbit. The best one 
can ask is that 2ir is the average of the distances around 
the various orbits. 

To consider this in detail, introduce a coordinate sys- 
tem, topologically the same as the standard spherical co- 
ordinates on the metric sphere, but adapted to the poten- 
tial function z so that the latitude lines are the level sur- 



11 In fact, using the higher eigenvectors, one could in principle com- 
pute higher-order multipole moments. We see this as a natural 
extension of the method laid out in [75] for defining the higher 
multipole moments of axisymmctric black holes. 



faces of z (and, in particular, the poles are at the two crit- 
ical points we have assumed z to have). More precisely, 
choose z for the zenith coordinate on the sphere, and an 
arbitrary rotational coordinate — say, the azimuthal an- 
gle in the encompassing spatial slice, describing rotations 
about the axis connecting the critical points of z — for the 
azimuthal coordinate ip on the sphere. If the parameter r 
is defined such that </> = (d/dT) z=cons t., then in the basis 
related to these coordinates, the components of 4> are: 

<P{z,<p) = (%) =0, (A14) 



if 



(A15) 



Around a closed orbit C(z), at fixed z, the parameter r 
changes by a value of: 



r(z) 



dp 



C{z) <t> v (z,tp) 
dip 
e<P z d z 

V§d(p, 



C{z) 



(A16) 
(A17) 
(A18) 



where g is the determinant of the surface metric, eval- 
uated in the (z, (p) coordinates. Note that Eq. (A18) 
follows from Eq. (A17) by the fact that the condition 
gABgcD£ AC £ BD = 2 implies e vz — \ j\fg. The average 
value of r, over the various orbits, is: 

sfldydz (A19) 
(A20) 



^max ^min J z ra i n J C{z) 

A 



^max ^min 



where A is the surface area of the apparent horizon. Re- 
quiring this average to equal 2w, wc arrive at the normal- 
ization condition: 



27f(£max ^min) A. 



(A21) 



This normalization condition requires finding the min- 
imum and maximum values of the function z, which is 
only computed on a discrete grid. In our spectral code, in 
particular, this numerical grid is quite coarse, so numeri- 
cal interpolation is needed, in combination with an opti- 
mization routine. We have implemented such routines to 
search for z m ; n and z max , but a numerically-cheaper nor- 
malization condition would be of interest. Such a condi- 
tion arises when one assumes that the black hole under 
consideration is approximately Kerr. In the Kerr metric, 
for the function z generating the true rotation generator 
of the Kerr horizon, the following identity holds: 



H 



(z-((z))fdA = ^ 



(A22) 



2G 



10" 



* 10 

a 

u 
o 

tl in" 9 



10 



12 



1 1 1 

: I^K 1 } 1 1 


□ 


i 

1 1 1 1 ' 


^^s. ^A 






— 65 ^ ^^-^ 












- □ ■ ■ □ X coord \^ 




v-vx sc 












L x — x X AKV ( norma li ze d by extrema) 




Q-Q X AKV (normalized by integral) 
i 


1 



10 



AH 



15 



FIG. 20: Error, relative to the analytic solution, of the spin 
on the horizon of a Kerr black hole in slightly deformed co- 
ordinates. The vertical axis represents |Xcom P utcd - Xanalyticl, 
and data are shown for the spin computed with the standard 
coordinate rotation vector (in deformed coordinates, so not 
a true Killing vector) , and with our approximate Killing vec- 
tors (AKV) using both the extremum norm, Eq. (A21), and 
the integral norm, Eq. (A22). The spin computed from the 
coordinate rotation vector quickly converges to a physically 
inaccurate result. The spin from approximate Killing vectors 
converges in resolution Lah to the correct value x — 1/2- 
Curves are also shown for the two spin measures defined in 
the Appendix B. These spin measures also converge exponen- 
tially to the physically correct result. 



where ((z)) is the average of z over the sphere. The exis- 
tence of an identity of this form is somewhat nontrivial: 
the fact that the right side is given purely by the hori- 
zon area, and that it does not involve the spin of the Kerr 
hole, is what makes this identity useful as a normalization 
condition. This normalization is much easier to impose, 
and requires significantly less numerical effort. 

To close the discussion of spin computed from approxi- 
mate Killing vectors, we demonstrate the effectiveness of 
the method in a simple test case: an analytic Kerr black 
hole in slightly deformed coordinates. We begin with a 
Kerr black hole of dimensionless spin parameter \ = 1/2, 
in Kerr-Schild coordinates, but we rescale the x-axis by a 
factor of 1.1. This rescaling of the x coordinate causes the 
coordinate rotation vector xd y — yd x to no longer be the 
true, geometrical rotation generator. And indeed, when 
we compute the quasilocal angular momentum (Al) on 
the horizon using this coordinate vector, the result con- 
verges to a physically inaccurate value, as demonstrated 
by the black dotted curve in Fig. 20. If, however, the ap- 
proximate Killing vectors described above are used, the 
result is not only convergent, but physically accurate. 
Because the accuracy is slightly better with the normal- 
ization condition of Eq. (A22), that is the condition we 
use for all results presented in this paper. 



APPENDIX B: SCALAR-CURVATURE SPIN (SC 
SPIN) 

In this appendix, we define a spin measure in terms of 
the intrinsic geometry of the horizon, which we compare 
with the AKV spin in Sec. V. The AKV spin described 
in Appendix A is a well-defined measure of black hole 
spin, even when the holes' horizons have only approxi- 
mate symmetries. At times sufficiently before or after 
the holes merge, however, the horizons will not be too 
tidally distorted and thus will not be too different from 
the exactly- axisymmetric horizons of Kerr black holes. 

By assuming that the geometric properties of the hori- 
zon behave precisely as they do for a Kerr black hole, 
one can infer the hole's spin from those properties. For 
instance, it is common to measure polar and equatorial 
circumferences of the apparent horizon; the spin is then 
obtained by finding the Kerr spacetime with the same 
circumferences [80-82] . 

To avoid introducing coordinate dependence by defin- 
ing "polar" and "equatorial" planes, we infer the spin 
from the horizon's intrinsic scalar curvature R. The hori- 
zon scalar curvature R has previously been studied ana- 
lytically for Kerr-Newman black holes [83] and for Kerr 
black holes perturbed by a distant moon [84] . Numerical 
studies of R have focused attention on the quasinormal 
ringing of single, perturbed, black holes [80] as well as on 
the shape of the individual and common event horizons 
in Misner data [85]. To our knowledge, the scalar curva- 
ture R has not been previously used to infer the horizon 
spin in numerical simulations. 

At a given point on a Kerr black hole's horizon, the 
horizon scalar curvature R depends only on the hole's 
mass M and spin S. The extrema of R can be expressed 
in terms of the irreducible mass and dimensionless spin 
of the Kerr black hole via Eqs. (l)-(2) as 



min(_R) 
max(_R) 



-1 + 2^1^ 



2Ml 
2 



(Bla) 



MU 4 



-2 + /Y 2 + 2 v / l-X 2 ) (Bib) 



Solving for x an d requiring it to be real yields x as a 
function of M- llr and either min(R) or max(i?). We take 
these functions as definitions of the spin, even when the 
space-time is not precisely Kerr: 



W) := 1 - 



M^ r min(i?) 



2 + 2^/2M i 2 rr max(i?) 
M£ r max( J R) 



(B2a) 



(B2b) 



The definitions of the spin given by Eqs. (B2a)-(B2b) 
are manifestly independent of spatial coordinates and 
are well-defined for black holes that are tidally deformed. 
Also, as they only involve the intrinsic two-dimensional 
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geometry of the apparent horizon, they are also mani- 
festly independent of boost gauge, in the sense described 
in the previous appendix. 

We expect Xscf anc ^ Xsc* *° ^ e reasonable measures 
only if tidal forces can be neglected. Tidal forces scale 
with the cube of the separation of the holes; for binary 
with holes of equal mass M and separation d. tidal cou- 
pling is negligible when max(i?) — min(-R) 3> M/d 3 . 

We find it convenient to compute R from i) the scalar 
curvature R associated with the three-dimensional metric 
gij of the spatial slice E, and ii) the outward-pointing 
unit- vector field s 4 that is normal to H. This can by done 
by means of Gauss's equation [e.g., Eq. (D.51) of Ref. [86] 
(note that the Riemann tensor in Ref. [86] disagrees with 
ours by an overall sign)] 

R = R — 2R ij sV - K 2 + K ij k i:j , (B3) 



where Rij and R were defined after Eq. (11), and where 
Kij denotes the extrinsic curvature of the the apparent 
horizon H embedded in E (not to be confused with , 
the extrinsic curvature of the slice E embedded in M.). 
The horizon extrinsic curvature is given by 

=ViSj- SiS k V kSj . (B4) 

Inserting Eq. (B4) into Eq. (B3) shows that R can be 
evaluated exclusively in terms of quantities defined on 
the three-dimensional spatial slice E. 

The accuracy of these spin measures is demonstrated 
in Fig. 20, which shows a Kerr black hole with x = 1/2 
in slightly deformed coordinates so that the coordinate 
rotation vector no longer generates a symmetry. Again, 
both Xscf an d Xsc* converge exponentially to the phys- 
ically accurate result. 
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